references: [1] E. U. Condon and H. Odabasi, "Atomic Structure". [2] T. Inui, Y. Tanabe, and Y. Onodera, "応用群論". [3] H. Kamimura, S, Sugano, and Y. Tanabe, "配位子場理論とその応用". (1) 定義 Off[General::spell1];Off[ClebschGordan::tri];Off[General::spell]; Off[ClebschGordan::phy]; (* 1.1 quantum number of the angular momentum of the 3d orbit. *) l3d=2; (* 1.2 parent {L, S, seniority} and |cfp|^2 {{L, S, seniority},{cfp,cfp,....}}. Table 5.2 of [1]. The sign of Sqrt[|cfp|^2] is indicated. *) par3d[1]={{0,0,0}}; cfp3d[1]={{{2,1/2,1},{1}} }; spingrp[1]={{1}}; par3d[2]={{2,1/2,1}}; cfp3d[2]={ {{0,0,0},{1}}, {{1,1,2},{1}}, {{2,0,2},{1}}, {{3,1,2},{1}}, {{4,0,2},{1}}}; spingrp[2]={{1,3,5},{2,4}}; par3d[3]={{0,0,0},{1,1,2},{2,0,2},{3,1,2},{4,0,2}}; cfp3d[3]={ {{1,1/2,3},{0,7,15,-8,0}/30}, {{1,3/2,3},{0,-8,0,-7,0}/15}, {{2,1/2,1},{16,-9,-5,-21,-9}/60}, {{2,1/2,3},{0,-49,45,21,-25}/140}, {{3,1/2,3},{0,28,-10,7,-25}/70}, {{3,3/2,3},{0,-1,0,4,0}/5}, {{4,1/2,3},{0,0,-10,21,11}/42}, {{5,1/2,3},{0,0,0,-1,1}/2} }; spingrp[3]={{1,3,4,5,7,8},{2,6}}; par3d[4]={{1,1/2,3},{1,3/2,3},{2,1/2,1},{2,1/2,3}, {3,1/2,3},{3,3/2,3},{4,1/2,3},{5,1/2,3}}; cfp3d[4]={ {{0,0,0},{0,0,1,0,0,0,0,0}/1}, {{0,0,4},{0,0,0,1,0,0,0,0}/1}, {{1,1,2},{-14,-64,135,-35,-56,-56,0,0}/360}, {{1,1,4},{25,-14,0,10,-25,16,0,0}/90}, {{2,0,2},{-42,0,105,45,28,0,-60,0}/280}, {{2,0,4},{42,0,0,20,63,0,15,0}/140}, {{2,1,4},{-14,49,0,60,-21,-21,45,0}/210}, {{2,2,4},{0,3,0,0,0,7,0,0}/10}, {{3,0,4},{120,0,0,200,-105,0,-3,-132}/560}, {{3,1,2},{16,-56,315,15,-14,224,90,110}/840}, {{3,1,4},{-200,-448,0,120,-175,-112,-405,220}/1680}, {{4,0,2},{0,0,189,-25,70,0,66,-154}/504}, {{4,0,4},{0,0,0,88,385,0,-507,-28}/1008}, {{4,1,4},{0,0,0,200,315,-560,297,308}/1680}, {{5,1,4},{0,0,0,0,5,20,-9,26}/60}, {{6,0,4},{0,0,0,0,0,0,3,7}/10} }; spingrp[4]={{1,2,5,6,9,12,13,16},{3,4,7,10,11,14,15},{8}}; par3d[5]={{0,0,0},{0,0,4},{1,1,2},{1,1,4},{2,0,2}, {2,0,4},{2,1,4},{2,2,4},{3,0,4},{3,1,2},{3,1,4}, {4,0,2},{4,0,4},{4,1,4},{5,1,4},{6,0,4}}; cfp3d[5]={ {{0,1/2,5},{0,0,0,0,0,-2,3,0,0,0,0,0,0,0,0,0}/5}, {{0,5/2,5},{0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0}/1}, {{1,1/2,3},{0,0,14,25,30,15,10,0,-15,-16,-25,0,0,0,0,0}/150}, {{1,3/2,3},{0,0,-64,14,0,0,35,-75,0,-56,56,0,0,0,0,0}/300}, {{2,1/2,1},{6,0,-9,0,-5,0,0,0,0,-21,0,-9,0,0,0,0}/50}, {{2,1/2,3},{0,-14,-49,-14,45,-10,60,0,35,21,-21,-25,-11,45,0,0}/350}, {{2,1/2,5},{0,-56,0,126,0,90,60,0,35,0,189,0,99,45,0,0}/700}, {{2,3/2,5},{0,0,0,126,0,0,-135,-175,0,0,-84,0,0,180,0,0}/700}, {{3,1/2,3},{0,0,448,-200,-160,180,120,0,105,112,-175,-400,275,-405,220,0}/2800}, {{3,1/2,5},{0,0,0,360,0,-100,600,0,-525,0,-315,0,495,-9,-396,0}/2800}, {{3,3/2,3},{0,0,-56,-16,0,0,-15,-175,0,224,14,0,0,-90,-110,0}/700}, {{4,1/2,3},{0,0,0,0,-800,-100,600,0,-7,1680,945,880,845,891,924,-728}/8400}, {{4,1/2,5},{0,0,0,0,0,1452,968,0,-2541,0,4235,0,-1215,-5577,-308,-2184}/18480}, {{4,3/2,5},{0,0,0,0,0,0,25,-105,0,0,-70,0,0,-66,154,0}/420}, {{5,1/2,3},{0,0,0,0,0,0,0,0,33,-220,55,220,-5,-99,286,182}/1100}, {{6,1/2,5},{0,0,0,0,0,0,0,0,0,0,0,0,-45,99,231,-175}/550} }; spingrp[5]={{1,3,5,6,7,9,10,12,13,15,16},{4,8,11,14},{2}}; (* 1.3 parent and cfp's for n>more than half. 5.14-(17) of [1], (6.178) of [2]. *) Module[{nfull,n,L,S,Ld,Sd},nfull=2*(2*l3d+1); Do[n=nfull-nn; par3d[nn]=Table[cfp3d[n+1][[i,1]],{i,1,Length[cfp3d[n+1]]}]; cfp3d[nn]=Table[ {L=par3d[n+1][[i,1]];S=par3d[n+1][[i,2]];par3d[n+1][[i]], Table[Ld=cfp3d[n+1][[j,1,1]];Sd=cfp3d[n+1][[j,1,2]]; (-1)^(L+S+Ld+Sd-l3d-1/2)*Sign[cfp3d[n+1][[j,2,i]]]* (n+1)*(2*Ld+1)*(2*Sd+1)/(4*l3d+2-n)/(2*L+1)/(2*S+1) *Abs[cfp3d[n+1][[j,2,i]]] ,{j,1,Length[cfp3d[n+1]]}]},{i,1,Length[par3d[n+1]]}]; spingrp[nn]=spingrp[n] ,{nn,2*l3d+2,nfull-1}] ]; (* 1.4 matrix element of . (6.111) of [2], (2.57) of [3], 3.7-(20') of [1]. *) ccc[l_,m_,ld_,md_,k_]:=((-1)^m*Sqrt[(2*l+1)*(2*ld+1)] *ThreeJSymbol[{l,0},{ld,0},{k,0}] *ThreeJSymbol[{l,-m},{ld,md},{k,m-md}]) (* 1.5 reduced matrix element . (6.110) of [2]. *) rmc[l_,ld_,k_]:=(-1)^l*Sqrt[(2*l+1)*(2*ld+1)]*ThreeJSymbol[{l,0},{ld,0},{k,0}]; (* 1.6 reduced matrix elements of U. 10.5-(8) of [1]. This is necessary to calculate the Coulomb pair interaction energy. *) rmtablU[l_,nel_,k_]:=Module[{Sa,Sb,La,Lb,L1}, Table[ Sa=cfp3d[nel][[i,1,2]]; La=cfp3d[nel][[i,1,1]]; Sb=cfp3d[nel][[j,1,2]]; Lb=cfp3d[nel][[j,1,1]]; KroneckerDelta[Sa,Sb]*nel*Sqrt[(2*La+1)*(2*Lb+1)] *Sum[L1=par3d[nel][[ii,1]];(-1)^(La+L1+l+k) *Sign[cfp3d[nel][[i,2,ii]]]*Sqrt[Abs[cfp3d[nel][[i,2,ii]]]] *Sign[cfp3d[nel][[j,2,ii]]]*Sqrt[Abs[cfp3d[nel][[j,2,ii]]]] *SixJSymbol[{La,Lb,k},{l,l,L1}] ,{ii,1,Length[par3d[nel]]}] ,{i,1,Length[cfp3d[nel]]},{j,1,Length[cfp3d[nel]]}]] (* 1.7 reduced matrix elements of V. 10.5-(10) of [1]. This is necessary to calculate the matrix elements of neutron scattering functions. *) rmtablV[l_,nel_,k_]:=Module[{Sa,Sb,La,Lb,L1,S1}, Table[ Sa=cfp3d[nel][[i,1,2]]; La=cfp3d[nel][[i,1,1]]; Sb=cfp3d[nel][[j,1,2]]; Lb=cfp3d[nel][[j,1,1]]; nel*Sqrt[(2*La+1)*(2*Lb+1)*(2*Sa+1)*(2*Sb+1)*3/2] *Sum[L1=par3d[nel][[ii,1]];S1=par3d[nel][[ii,2]];(-1)^(La+L1+l+k+Sa+S1+3/2) *Sign[cfp3d[nel][[i,2,ii]]]*Sqrt[Abs[cfp3d[nel][[i,2,ii]]]] *Sign[cfp3d[nel][[j,2,ii]]]*Sqrt[Abs[cfp3d[nel][[j,2,ii]]]] *SixJSymbol[{Sa,Sb,1},{1/2,1/2,S1}]*SixJSymbol[{La,Lb,k},{l,l,L1}] ,{ii,1,Length[par3d[nel]]}] ,{i,1,Length[cfp3d[nel]]},{j,1,Length[cfp3d[nel]]}]] (* 1.8 A routine to check if the obtained eigenvectors are diagonalized. *) digcheck[vecs_]:=Table[Chop[N[ Sum[Conjugate[vecs[[i,k]]]*vecs[[j,k]],{k,1,Length[vecs]}]]] ,{i,1,Length[vecs]},{j,1,Length[vecs]}] (* 1.9 Definition of the raising and lowering operators. *) jp[j_,m_,n_]:=(Product[Sqrt[(j-k) (j+k+1)],{k,m,m+n-1,1}]); jm[j_,m_,n_]:=(Product[Sqrt[(j+k) (j-k+1)],{k,m,m-n+1,-1}]); (2) d電子数とクーロン相互作用パラメータ（Slater積分） (* 2.1 Input the number of 3d electrons and Slater integrals. *) n3d=2; f2=10; f4=6.25; rmV11=rmtablV[l3d,n3d,1]; (* 2.2 matrix elements of Coulomb interaction energy. 5.10-(19) of [1]. F[0], F[2], and F[4] are the Slater integrals. *) utabl=Table[rmtablU[l3d,n3d,kk],{kk,0,2*l3d,2}]; coultabl=Module[{La,Sa,Lb,Sb,L2},Table[ Sum[ La=cfp3d[n3d][[i,1,1]];Sa=cfp3d[n3d][[i,1,2]]; Lb=cfp3d[n3d][[j,1,1]];Sb=cfp3d[n3d][[j,1,2]]; F[k]*Chop[N[ Abs[rmc[l3d,l3d,k]]^2/2 *(1/(2*La+1)*Sum[(-1)^(La+cfp3d[n3d][[n,1,1]]) *utabl[[k/2+1,i,n]]*utabl[[k/2+1,n,j]] *KroneckerDelta[Sa,Sb]*KroneckerDelta[La,Lb] ,{n,1,Length[cfp3d[n3d]]}] -n3d/(2*l3d+1)*KroneckerDelta[i,j]) ]],{k,0,2*l3d,2}] ,{i,1,Length[cfp3d[n3d]]},{j,1,Length[cfp3d[n3d]]}]]; Print["Coulomb interaction matrix element: "] Print[Tab,coultabl] (* 2.3 Configuration average of the Coulomb pair interaction energy. See 5.10-(22) of [1]. *) ave3d=Expand[ N[n3d*(n3d-1)/2]*(F[0]-N[2/63]*F[2]-N[2/63]*F[4]) ]; Print["Average Coulomb interaction energy: "] Print[Tab,ave3d] (* 2.4 Matrix elements of the Coulomb interaction energy relative to the configuration average. See 5.10-(22) and appendix 3 of [1].*) couldiff=Chop[Table[If[i==j,coultabl[[i,j]]-ave3d,coultabl[[i,j]]], {i,1,Length[cfp3d[n3d]]},{j,1,Length[cfp3d[n3d]]}]]; Print["Difference from average Coulomb interaction: "] Print[Tab,couldiff] (* 2.5 Matrix elements of the Coulomb pair interaction. Configuration average is removed from diagonal elements. *) coulelem=Module[{L,Ld}, Table[ L=cfp3d[n3d][[i,1,1]];Ld=cfp3d[n3d][[j,1,1]]; Table[ If[ML==MLd,couldiff[[i,j]]/.{F[2]->f2,F[4]->f4},0], {ML,L,-L,-1},{MLd,Ld,-Ld,-1} ], {i,1,Length[cfp3d[n3d]]},{j,1,Length[cfp3d[n3d]]}] ]; Print["Average Coulomb interaction (eV): "] Print[Tab,MatrixForm[ave3d] /.{F[2]->f2,F[4]->f4}] Print["Difference from average Coulomb interaction (eV): "] Print[Tab,MatrixForm[couldiff] /.{F[2]->f2,F[4]->f4}] (3) 点電荷モデルによる結晶場ハミルトニアン (* 3.1 Input anion sites, charges, and the averages of the radial expansions. *) pcsites={{0,0,1},{0,0,-1},{1,0,0},{-1,0,0},{0,1,0},{0,-1,0}}; charge={-1,-1,-1,-1,-1,-1}*1.5; rave[2]=0.446; rave[4]=0.457; (* 3.2 calculation of the CEF parameters. *) 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=N[ArcTan[xd,yd]]]; Which[zd>0 && xd==0 && yd==0,theta=0,zd<0 && xd==0 && yd==0,theta=Pi,True, theta=N[ArcCos[zd/r]]]; r,theta,phi},{i,1,Length[pcsites]}]]; Print["Atomic Positions (r,theta,phi)="]; Print[pctp]; Do[Do[qfac[k,m]=Sqrt[4*Pi/(2*k+1)]*Chop[Sum[charge[[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,4}]; Print["Expansion coefficient of the CEF potential. q[odd,m] becomes 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["When octahedral cubic symmetry: 10Dq=",qfac[4,0]*2/7*10*rave[4],"(eV)"]; (* 3.3 reduced matrix elements for CEF operators. (8.124) of [3]. *) rmcef[0]=rmc[l3d,l3d,0]*utabl[[1]]; rmcef[2]=rmc[l3d,l3d,2]*utabl[[2]]; rmcef[4]=rmc[l3d,l3d,4]*utabl[[3]]; (* 3.4 Matrix elements of the CEF Hamiltonian. (8.48) of [3]. *) cefelem=Module[{L,Ld}, Table[ L=cfp3d[n3d][[i,1,1]];Ld=cfp3d[n3d][[j,1,1]]; Table[ Sum[Sum[ If[q==(ML-MLd) && qfac[k,q]!=0 && Abs[L-Ld]<=k<=(L+Ld), (-rave[k]*6*qfac[k,q]*rmcef[k][[i,j]]* ClebschGordan[{Ld,MLd},{k,q},{L,ML}]/Sqrt[2*L+1]),0 ] ,{q,-k,k}],{k,2,4,2}] ,{ML,L,-L,-1},{MLd,Ld,-Ld,-1}], {i,1,Length[cfp3d[n3d]]},{j,1,Length[cfp3d[n3d]]}] ]; (* 3.5 Matrix elements of very weak zeeman energy. This is necessary to obtain orthogonal eigenvectors. *) weakz=10^(-8); zeemelem=Module[{L,Ld,S,Sd}, Table[L=cfp3d[n3d][[i,1,1]];Ld=cfp3d[n3d][[j,1,1]]; S=cfp3d[n3d][[i,1,2]];Sd=cfp3d[n3d][[j,1,2]]; Table[ If[L==Ld && S==Sd && ML==MLd,ML*weakz,0], {ML,L,-L,-1},{MLd,Ld,-Ld,-1} ], {i,1,Length[cfp3d[n3d]]},{j,1,Length[cfp3d[n3d]]}] ]; (* 3.6 Show groups of multiplets with same spin quantum number *) Do[Print[i,Tab,"S=",cfp3d[n3d][[spingrp[n3d][[i,1]],1,2]],Tab, "Multiplet Group=",spingrp[n3d][[i]] ] ,{i,1,Length[spingrp[n3d]]}] (4) クーロン相互作用＋結晶場ハミルトニアンの固有関数 (* 4.1 Select one of the multiplet groups. The CEF eigenstates are calculated. *) mnum=2; matgrp=spingrp[n3d][[mnum]]; dcef=N[Flatten[Table[Table[ Flatten[Table[cefelem[[matgrp[[i]],matgrp[[j]]]][[k]],{j,1,Length[matgrp]}]] ,{k,1,Length[cefelem[[matgrp[[i]],matgrp[[1]]]]]}] ,{i,1,Length[matgrp]}],1] ]; dcoul=N[Flatten[Table[Table[ Flatten[Table[coulelem[[matgrp[[i]],matgrp[[j]]]][[k]],{j,1,Length[matgrp]}]] ,{k,1,Length[coulelem[[matgrp[[i]],matgrp[[1]]]]]}],{i,1,Length[matgrp]}],1] ]; dzeem=N[Flatten[Table[Table[ Flatten[Table[zeemelem[[matgrp[[i]],matgrp[[j]]]][[k]],{j,1,Length[matgrp]}]] ,{k,1,Length[zeemelem[[matgrp[[i]],matgrp[[1]]]]]}],{i,1,Length[matgrp]}],1] ]; {vals,vecs}=Eigensystem[N[dcef+dcoul+dzeem]]; {vals,vecs}=Chop[{vals,vecs}]; valsvecs=Table[{vals[[i]],vecs[[i]]},{i,1,Length[vals]}]; valsvecs=Sort[valsvecs]; vals=Table[Chop[valsvecs[[i,1]]],{i,1,Length[vals]}]; vecs=Table[Chop[valsvecs[[i,2]]],{i,1,Length[vals]}]; cfvecs=Module[{i,L},Table[i=0;Flatten[Table[ L=cfp3d[n3d][[matgrp[[ii]],1,1]]; Table[i=i+1;{vecs[[hh,i]],{L,L-j},matgrp[[ii]]},{j,0,2*L}], {ii,1,Length[matgrp]}],1],{hh,1,Length[vecs]}]]; cfvals=vals; cfvecs2=Module[{i,dummy},Table[dummy={};i=1; While[i<=Length[cfvecs[[hh]]], If[Chop[cfvecs[[hh,i,1]],10^(-6)]==0,0,dummy= Join[dummy,{cfvecs[[hh,i]]}]]; i++];dummy,{hh,1,Length[cfvecs]}]]; Print["S=",cfp3d[n3d][[matgrp[[1]],1,2]],Tab, "L=",Table[cfp3d[n3d][[matgrp[[i]],1,1]],{i,1,Length[matgrp]}]]; Print["Eigenvalues (eV) and eigenvectors of CEF+Coulomb Hamiltonian:"] Do[Print[i,Tab,cfvals[[i]],Chop[cfvecs2[[i]],10^(-6)]],{i,1,Length[cfvals]}] (* 4.2 Check if the obtained eigenvectors are orthogonalized. *) Chop[ digcheck[vecs], weakz*4] Table[Chop[(dcef + dcoul + dzeem).vecs[[i]] - vals[[i]]*vecs[[i]], weakz*4], {i, 1, Length[vals]}] (5) スピン軌道相互作用とゼーマン効果 (* 5.1 Select the Coulomb+CEF eigen states. The matrix elements of L and S are calculated. *) ii=1; if=3; gsvals=Table[cfvals[[i]],{i,ii,if}]; gsvecs=Table[cfvecs[[i]],{i,ii,if}]; gsvecs2=Module[{i,dummy},Table[dummy={};i=1; While[i<=Length[gsvecs[[hh]]], If[Chop[gsvecs[[hh,i,1]],10^(-6)]==0,0,dummy= Join[dummy,{gsvecs[[hh,i]]}]]; i++];dummy,{hh,1,Length[gsvecs]}]]; Print["Selected CEF+Coulomb eigenenergies (eV) and eigenvectors: "] Do[Print[gsvals[[i]],Chop[gsvecs2[[i]],10^(-6)]],{i,1,Length[gsvals]}] Clear[ii,if] Module[{S,L,MS,ML,Ld,MLd,fac}, S=cfp3d[n3d][[matgrp[[1]],1,2]]; lxelem=Flatten[Table[Table[Flatten[Table[ Sum[L=gsvecs[[i,ii,2,1]];ML=gsvecs[[i,ii,2,2]]; Ld=gsvecs[[j,jj,2,1]];MLd=gsvecs[[j,jj,2,2]]; fac=Conjugate[gsvecs[[i,ii,1]]]*gsvecs[[j,jj,1]]; If[fac!=0 && L==Ld, fac*(If[MS==MSd,Switch[ML-MLd,1,jp[L,MLd,1]/2,-1,jm[L,MLd,1]/2,_,0],0]) , 0], {ii,1,Length[gsvecs[[1]]]},{jj,1,Length[gsvecs[[1]]]} ] ,{MSd,S,-S,-1},{j,1,Length[gsvecs]}]],{i,1,Length[gsvecs]}],{MS,S,-S,-1}],1]; sxelem=Flatten[Table[Table[Flatten[Table[ Sum[L=gsvecs[[i,ii,2,1]];ML=gsvecs[[i,ii,2,2]]; Ld=gsvecs[[j,jj,2,1]];MLd=gsvecs[[j,jj,2,2]]; fac=Chop[Conjugate[gsvecs[[i,ii,1]]]*gsvecs[[j,jj,1]]]; If[fac!=0 && L==Ld, fac*(If[ML==MLd,Switch[MS-MSd,1,jp[S,MSd,1]/2,-1,jm[S,MSd,1]/2,_,0],0]) , 0], {ii,1,Length[gsvecs[[1]]]},{jj,1,Length[gsvecs[[1]]]} ] ,{MSd,S,-S,-1},{j,1,Length[gsvecs]}]],{i,1,Length[gsvecs]}],{MS,S,-S,-1}],1]; lyelem=Flatten[Table[Table[Flatten[Table[ Sum[L=gsvecs[[i,ii,2,1]];ML=gsvecs[[i,ii,2,2]]; Ld=gsvecs[[j,jj,2,1]];MLd=gsvecs[[j,jj,2,2]]; fac=Conjugate[gsvecs[[i,ii,1]]]*gsvecs[[j,jj,1]]; If[fac!=0 && L==Ld, fac*(If[MS==MSd,Switch[ML-MLd,1,-I*jp[L,MLd,1]/2,-1,I*jm[L,MLd,1]/2,_,0],0]) , 0], {ii,1,Length[gsvecs[[1]]]},{jj,1,Length[gsvecs[[1]]]} ] ,{MSd,S,-S,-1},{j,1,Length[gsvecs]}]],{i,1,Length[gsvecs]}],{MS,S,-S,-1}],1]; syelem=Flatten[Table[Table[Flatten[Table[ Sum[L=gsvecs[[i,ii,2,1]];ML=gsvecs[[i,ii,2,2]]; Ld=gsvecs[[j,jj,2,1]];MLd=gsvecs[[j,jj,2,2]]; fac=Conjugate[gsvecs[[i,ii,1]]]*gsvecs[[j,jj,1]]; If[fac!=0 && L==Ld, fac*(If[ML==MLd,Switch[MS-MSd,1,-I*jp[S,MSd,1]/2,-1,I*jm[S,MSd,1]/2,_,0],0]) , 0], {ii,1,Length[gsvecs[[1]]]},{jj,1,Length[gsvecs[[1]]]} ] ,{MSd,S,-S,-1},{j,1,Length[gsvecs]}]],{i,1,Length[gsvecs]}],{MS,S,-S,-1}],1]; lzelem=Flatten[Table[Table[Flatten[Table[ Sum[L=gsvecs[[i,ii,2,1]];ML=gsvecs[[i,ii,2,2]]; Ld=gsvecs[[j,jj,2,1]];MLd=gsvecs[[j,jj,2,2]]; fac=Conjugate[gsvecs[[i,ii,1]]]*gsvecs[[j,jj,1]]; If[fac!=0 && L==Ld,fac*(If[MS==MSd && ML==MLd,ML,0]) , 0], {ii,1,Length[gsvecs[[1]]]},{jj,1,Length[gsvecs[[1]]]} ] ,{MSd,S,-S,-1},{j,1,Length[gsvecs]}]],{i,1,Length[gsvecs]}],{MS,S,-S,-1}],1]; szelem=Flatten[Table[Table[Flatten[Table[ Sum[L=gsvecs[[i,ii,2,1]];ML=gsvecs[[i,ii,2,2]]; Ld=gsvecs[[j,jj,2,1]];MLd=gsvecs[[j,jj,2,2]]; fac=Conjugate[gsvecs[[i,ii,1]]]*gsvecs[[j,jj,1]]; If[fac!=0 && L==Ld,fac*(If[ML==MLd && MS==MSd,MS,0]) , 0], {ii,1,Length[gsvecs[[1]]]},{jj,1,Length[gsvecs[[1]]]} ] ,{MSd,S,-S,-1},{j,1,Length[gsvecs]}]],{i,1,Length[gsvecs]}],{MS,S,-S,-1}],1]; gscfel=Flatten[Table[Table[Flatten[Table[ If[i==j && MS==MSd,(gsvals[[i]]-gsvals[[1]])*1000,0] ,{MSd,S,-S,-1},{j,1,Length[gsvecs]}]],{i,1,Length[gsvecs]}],{MS,S,-S,-1}],1]; ] (* 5.2 Set the spin-orbit coupling constant zetaso [meV], magnitude of the magnetic field hhh [T], and the unit vector hhhdirec of the direction of the field. Then the eigenstates are calculated. *) zetaso=10; hhh=0; hhhdirec={0,0,1}; magf=0.05796*hhh*hhhdirec; rmso=N[ zetaso*Sqrt[2*3*5]*rmV11 ]; Module[{S,L,MS,ML,Ld,MLd,q,fac}, S=cfp3d[n3d][[matgrp[[1]],1,2]]; soelem=Chop[Flatten[Table[Table[Flatten[Table[ Sum[L=gsvecs[[i,ii,2,1]];ML=gsvecs[[i,ii,2,2]]; Ld=gsvecs[[j,jj,2,1]];MLd=gsvecs[[j,jj,2,2]];q=(-ML+MLd); fac=Chop[(rmso[[gsvecs[[i,ii,3]],gsvecs[[j,jj,3]]]] *Conjugate[gsvecs[[i,ii,1]]]*gsvecs[[j,jj,1]])]; If[fac!=0 && q==(MS-MSd) && Abs[q]<=1, N[fac/Sqrt[2*L+1]/Sqrt[2*S+1]*(-1)^q*ClebschGordan[{S,MSd},{1,q},{S,MS}] *ClebschGordan[{Ld,MLd},{1,-q},{L,ML}] ] , 0], {ii,1,Length[gsvecs[[1]]]},{jj,1,Length[gsvecs[[1]]]} ] ,{MSd,S,-S,-1},{j,1,Length[gsvecs]}]],{i,1,Length[gsvecs]}],{MS,S,-S,-1}],1]]; ] zeeman=Chop[(lxelem+2*sxelem)*magf[[1]]+(lyelem+2*syelem)*magf[[2]] +(lzelem+2*szelem)*magf[[3]]]; Print["Spin-Orbit=",MatrixForm[soelem]] Print["Zeeman=",MatrixForm[zeeman]] {vals2,vecs2}=Chop[Eigensystem[gscfel+soelem+zeeman]]; valsvecs=Table[{vals2[[i]],vecs2[[i]]},{i,1,Length[vals2]}]; valsvecs=Sort[valsvecs]; vals2=Table[Chop[valsvecs[[i,1]]],{i,1,Length[vals2]}]; vecs2=Table[Chop[valsvecs[[i,2]]],{i,1,Length[vals2]}]; Print["Eigen energies (meV) and eigenvectors: "] Do[Print[vals2[[i]],Chop[vecs2[[i]],10^(-6)]],{i,1,Length[vecs2]}] Module[{L,ML,S,MS,mulnum,i,dummy}, S=cfp3d[n3d][[matgrp[[1]],1,2]]; sovecs=Table[Flatten[Table[iii=0; L=gsvecs[[1,j,2,1]];ML=gsvecs[[1,j,2,2]];mulnum=gsvecs[[1,j,3]]; Table[{Sum[iii=iii+1;vecs2[[ii,iii]]*gsvecs[[i,j,1]],{i,1,Length[gsvecs]}] ,{L,ML,S,MS},mulnum},{MS,S,-S,-1}] ,{j,1,Length[gsvecs[[1]]]}],1],{ii,1,Length[vecs2]}]; sovecs2=Table[dummy={};i=1; While[i<=Length[sovecs[[hh]]], If[Chop[sovecs[[hh,i,1]],10^(-6)]==0,0,dummy= Join[dummy,{sovecs[[hh,i]]}]]; i++];dummy,{hh,1,Length[sovecs]}];] sovals=vals2; Print[] Print["Eigen States in {L,ML,S,MS} representation: "] Do[Print[sovals[[i]],Chop[sovecs2[[i]],10^(-6)]],{i,1,Length[sovecs2]}] (* 5.3 Set the temperature in Kelvin. The magnetization is calculated. *) temp=0.1; za=Sum[Exp[-sovals[[i]]*11.6/temp],{i,1,Length[sovals]}]; mag=Module[{L,ML,Ld,MLd,S,MS,Sd,MSd}, Sum[Sum[L=sovecs2[[i,ii,2,1]] ;ML=sovecs2[[i,ii,2,2]]; S=sovecs2[[i,ii,2,3]] ;MS=sovecs2[[i,ii,2,4]]; Ld=sovecs2[[i,jj,2,1]] ;MLd=sovecs2[[i,jj,2,2]]; Sd=sovecs2[[i,jj,2,3]] ;MSd=sovecs2[[i,jj,2,4]]; fac=Conjugate[sovecs2[[i,ii,1]]]*sovecs2[[i,jj,1]]; -{fac*(If[MS==MSd && L==Ld,Switch[ML-MLd,1,jp[L,MLd,1]/2,-1,jm[L,MLd,1]/2,_,0],0] +2*If[ML==MLd && L==Ld,Switch[MS-MSd,1,jp[S,MSd,1]/2,-1,jm[S,MSd,1]/2,_,0],0]), fac*(If[MS==MSd && L==Ld,Switch[ML-MLd,1,-I*jp[L,MLd,1]/2,-1,I*jm[L,MLd,1]/2,_,0],0] +2*If[ML==MLd && L==Ld,Switch[MS-MSd,1,-I*jp[S,MSd,1]/2,-1,I*jm[S,MSd,1]/2,_,0],0]), fac*If[MS==MSd && ML==MLd && L==Ld,ML+2*MS,0]} ,{ii,1,Length[sovecs2[[i]]]},{jj,1,Length[sovecs2[[i]]]}] *Exp[-sovals[[i]]*11.6/temp] ,{i,1,Length[sovals]}]/za]; Print["Magnetization=",Chop[mag]," (muB)"] (* 5.4 Matrix elements of L and S for the eigenstates. *) Module[{L,ML,Ld,MLd,S,MS,Sd,MSd},elem2=Table[ Sum[L=sovecs2[[i,ii,2,1]] ;ML=sovecs2[[i,ii,2,2]]; S=sovecs2[[i,ii,2,3]] ;MS=sovecs2[[i,ii,2,4]]; Ld=sovecs2[[j,jj,2,1]] ;MLd=sovecs2[[j,jj,2,2]]; Sd=sovecs2[[j,jj,2,3]] ;MSd=sovecs2[[j,jj,2,4]]; fac=Conjugate[sovecs2[[i,ii,1]]]*sovecs2[[j,jj,1]]; {fac*If[MS==MSd && L==Ld,Switch[ML-MLd,1,jp[L,MLd,1]/2,-1,jm[L,MLd,1]/2,_,0],0], fac*If[ML==MLd && L==Ld,Switch[MS-MSd,1,jp[S,MSd,1]/2,-1,jm[S,MSd,1]/2,_,0],0], fac*If[MS==MSd && L==Ld,Switch[ML-MLd,1,-I*jp[L,MLd,1]/2,-1,I*jm[L,MLd,1]/2,_,0],0], fac*If[ML==MLd && L==Ld,Switch[MS-MSd,1,-I*jp[S,MSd,1]/2,-1,I*jm[S,MSd,1]/2,_,0],0], fac*If[MS==MSd && ML==MLd && L==Ld,ML,0], fac*If[MS==MSd && ML==MLd && L==Ld,MS,0]} ,{ii,1,Length[sovecs2[[i]]]},{jj,1,Length[sovecs2[[j]]]}] ,{i,1,Length[sovecs2]},{j,1,Length[sovecs2]}]]; lxelem2=Table[elem2[[i,j,1]],{i,1,Length[sovecs2]},{j,1,Length[sovecs2]}]; sxelem2=Table[elem2[[i,j,2]],{i,1,Length[sovecs2]},{j,1,Length[sovecs2]}]; lyelem2=Table[elem2[[i,j,3]],{i,1,Length[sovecs2]},{j,1,Length[sovecs2]}]; syelem2=Table[elem2[[i,j,4]],{i,1,Length[sovecs2]},{j,1,Length[sovecs2]}]; lzelem2=Table[elem2[[i,j,5]],{i,1,Length[sovecs2]},{j,1,Length[sovecs2]}]; szelem2=Table[elem2[[i,j,6]],{i,1,Length[sovecs2]},{j,1,Length[sovecs2]}]; (6) 磁気形状因子と中性子散乱 Reference: "Theory of Magnetic Neutron and Photon Scattering", E. Balcar and S. W. Lovesey, (Oxford, 1989) (* 6.1.1 Reduced matrix elements of U and V *) utablns=Table[Chop[rmtablU[l3d,n3d,kk]],{kk,1,2*l3d-1,2}]; vtablns=Table[Chop[rmtablV[l3d,n3d,kk]],{kk,0,2*l3d,2}]; (* 6.1.2 Definition (3.4.9) *) AA[K_,Kd_,l_]:=((-1)^(l+1)*Sqrt[(2*l+3)*(l+1)/(2*l+1)] *ThreeJSymbol[{l,0},{K,0},{l+1,0}]*SixJSymbol[{1,Kd,K},{l,l+1,l}]) (* 6.1.3 Matrix elements of orbital and spin part of the scattering operator Q. Eq. (3.9.3) and (3.9.9) *) orbtQ[v_,L_,ML_,S_,MS_,vd_,Ld_,MLd_,Sd_,MSd_,q_,th_,ph_]:= Sum[ If[ OddQ[k],-(rj[k-1]+rj[k+1])* N[Sqrt[4*Pi]*I^(k+1)*(2*k+1)/Sqrt[3*(k+1)]* (-1)^(ML+Ld-k)*(2*l3d+1)^2*AA[k,k,l3d]* ThreeJSymbol[{L,-ML},{Ld,MLd},{k,ML-MLd}]* SphericalHarmonicY[k-1,q-ML+MLd,th,ph]* ClebschGordan[{k-1,q-ML+MLd},{k,ML-MLd},{1,q}]* utablns[[(k+1)/2,v,vd]]],0] ,{k,Max[Abs[L-Ld],1],Min[L+Ld,2*l3d-1]}] spinQ[v_,L_,ML_,S_,MS_,vd_,Ld_,MLd_,Sd_,MSd_,q_,th_,ph_]:=If[q==MS-MSd, Sum[ If[EvenQ[k], rj[k]*N[Sqrt[4*Pi]*I^k*(2*l3d+1)*Sqrt[2*k+1]*(-1)^(MS-ML+L+l3d+S)* ThreeJSymbol[{l3d,0},{k,0},{l3d,0}]*vtablns[[k/2+1,v,vd]]* ThreeJSymbol[{Sd,MSd},{S,-MS},{1,MS-MSd}]* ThreeJSymbol[{Ld,MLd},{L,-ML},{k,ML-MLd}]* SphericalHarmonicY[k,-(ML-MLd),th,ph]*(-1)^(ML-MLd)],0], {k,Max[0,Abs[L-Ld]],Min[2*l3d,L+Ld]}],0] (* 6.1.4 Spin-Orbit and Zeeman eigen states obtained in sec. 5. *) vals=sovals;vecs=sovecs2; Do[Print[i,Tab,vals[[i]],Chop[vecs[[i]]]],{i,1,Length[vecs]}] (* Approximtion of , and . Input coefficients of the ion you want to calculate. *) rjcoef[0]={0.4198,14.283,0.6054,5.469,0.9241,-0.009,-0.9498}; rjcoef[2]={1.2427,14.997,1.9567,6.118,0.5732,2.258,0.0031}; rjcoef[4]={-0.6603,13.607,0.2322,6.218,0.4104,1.74,0.0101}; frj[l_,s_]:=((Sum[rjcoef[l][[i]]*Exp[-rjcoef[l][[i+1]]*s^2],{i,1,5,2}] +rjcoef[l][[7]])*If[l==0,1,s^2]); Plot[{frj[0,s],frj[2,s],frj[4,s]},{s,0,1.2},PlotRange->All ,AxesLabel->{"",""} ,PlotStyle->{RGBColor[0,0,0],RGBColor[1,0,0],RGBColor[0,0,1]}] (* 6.3.1 Select initial states and final states. *) iina=1;iinb=1; ifna=1;ifnb=1; (* 6.3.2 Table of the matrix elements of Q between initial and final states. *) Clear[rj0,rj2,rj4,th,ph]; Module[{v,L,ML,S,MS,vd,Ld,MLd,Sd,MSd}, Qelspin=Table[Table[ elem=Chop[Sum[ v=vecs[[i,ii,3]];L=vecs[[i,ii,2,1]];ML=vecs[[i,ii,2,2]]; S=vecs[[i,ii,2,3]];MS=vecs[[i,ii,2,4]]; vd=vecs[[j,jj,3]];Ld=vecs[[j,jj,2,1]];MLd=vecs[[j,jj,2,2]]; Sd=vecs[[j,jj,2,3]];MSd=vecs[[j,jj,2,4]]; Expand[Conjugate[vecs[[i,ii,1]]]*vecs[[j,jj,1]]* spinQ[v,L,ML,S,MS,vd,Ld,MLd,Sd,MSd,q,th,ph]] ,{ii,1,Length[vecs[[i]]]},{jj,1,Length[vecs[[j]]]}] ]; {Coefficient[elem,rj[0]],Coefficient[elem,rj[2]],Coefficient[elem,rj[4]]} ,{q,-1,1}] ,{i,iina,iinb},{j,ifna,ifnb}]; Qelorbt=Table[Table[ elem=Chop[Sum[ v=vecs[[i,ii,3]];L=vecs[[i,ii,2,1]];ML=vecs[[i,ii,2,2]]; S=vecs[[i,ii,2,3]];MS=vecs[[i,ii,2,4]]; vd=vecs[[j,jj,3]];Ld=vecs[[j,jj,2,1]];MLd=vecs[[j,jj,2,2]]; Sd=vecs[[j,jj,2,3]];MSd=vecs[[j,jj,2,4]]; Expand[Conjugate[vecs[[i,ii,1]]]*vecs[[j,jj,1]]* orbtQ[v,L,ML,S,MS,vd,Ld,MLd,Sd,MSd,q,th,ph]] ,{ii,1,Length[vecs[[i]]]},{jj,1,Length[vecs[[j]]]}] ]; {Coefficient[elem,rj[0]],Coefficient[elem,rj[2]],Coefficient[elem,rj[4]]} ,{q,-1,1}] ,{i,iina,iinb},{j,ifna,ifnb}]; ] (* 6.3.3 Transform Q from spherical to Cartesian components. *) Qxyzspin=Table[{(Qelspin[[i,j,1]]-Qelspin[[i,j,3]])/Sqrt[2], (Qelspin[[i,j,1]]+Qelspin[[i,j,3]])*I/Sqrt[2], Qelspin[[i,j,2]]},{i,1,iinb-iina+1},{j,1,ifnb-ifna+1}]; Qxyzorbt=Table[{(Qelorbt[[i,j,1]]-Qelorbt[[i,j,3]])/Sqrt[2], (Qelorbt[[i,j,1]]+Qelorbt[[i,j,3]])*I/Sqrt[2], Qelorbt[[i,j,2]]},{i,1,iinb-iina+1},{j,1,ifnb-ifna+1}]; (* 6.3.4 Function of Qxyz table between initial and final states. *) Qtablspin[th_,ph_,rj0_,rj2_,rj4_]:=Evaluate[Table[Table[Chop[ ComplexExpand[ Qxyzspin[[i,j,alp,1]]*rj0+Qxyzspin[[i,j,alp,2]]*rj2 +Qxyzspin[[i,j,alp,3]]*rj4] ],{alp,1,3}], {i,1,iinb-iina+1},{j,1,ifnb-ifna+1}]]; Qtablorbt[th_,ph_,rj0_,rj2_,rj4_]:=Evaluate[Table[Table[Chop[ ComplexExpand[ Qxyzorbt[[i,j,alp,1]]*rj0+Qxyzorbt[[i,j,alp,2]]*rj2 +Qxyzorbt[[i,j,alp,3]]*rj4] ],{alp,1,3}], {i,1,iinb-iina+1},{j,1,ifnb-ifna+1}]]; (* 6.3.5 Input scattering vector in 1/Angstrom unit. *) qvec={0.00001,0,0}; (* 6.3.6 Calculation of Q, Q_perp and S(q,w). *) qx=qvec[[1]];qy=qvec[[2]];qz=qvec[[3]]; tau=Sqrt[qvec[[1]]^2+qvec[[2]]^2+qvec[[3]]^2]; kapp=qvec/tau; qval=tau/4/Pi; Which[qx==0 && qy>0,phi=Pi/2,qx==0 && qy<0,phi=3*Pi/2,qy==0 && qx>0,phi=0, qy==0 && qx<0,phi=Pi,qx==0 && qy==0,phi=0,True,phi=N[ArcTan[qx,qy]]]; Which[qz>0 && qx==0 && qy==0,theta=0,qz<0 && qx==0 && qy==0,theta=Pi,True,theta=N[ArcCos[qz/tau]]]; fftablspin=Chop[N[ Qtablspin[theta,phi,frj[0,qval],frj[2,qval],frj[4,qval]] ]]; fftablorbt=Chop[N[ Qtablorbt[theta,phi,frj[0,qval],frj[2,qval],frj[4,qval]] ]]; fftabl=fftablspin+fftablorbt; qperp=Table[Cross[kapp,Cross[fftabl[[i,j]],kapp]] ,{i,1,iinb-iina+1},{j,1,ifnb-ifna+1}]; sqwtabl=Table[Chop[Sum[ Conjugate[qperp[[i,j,alp]]]*qperp[[i,j,alp]],{alp,1,3}]] , {i,1,iinb-iina+1},{j,1,ifnb-ifna+1}]; Print["Kappa/4Pi=",N[qval]," (1/Angstrom)",Tab, "theta=",theta*180/Pi,"(deg)",Tab,"phi=",phi*180/Pi,"(deg)"] Print["{Qx,Qy,Qz}_S=",MatrixForm[fftablspin]] Print["{Qx,Qy,Qz}_O=",MatrixForm[fftablorbt]] Print["{Qx,Qy,Qz}=",MatrixForm[fftabl]] Print["{Qx,Qy,Qz}_perp=",MatrixForm[qperp]] Print["for unpolarized incident neutrons:"] Print["S(Kappa,Omega)=",MatrixForm[sqwtabl]] (* 6.3.7 Temperature factor. *) temp=300; z=Sum[Exp[-vals[[i]]*11.6056/temp],{i,1,Length[vals]}]; tempfac=Table[Exp[-vals[[i]]*11.6056/temp]/z,{i,iina,iinb}]; sqwtabl2=Table[sqwtabl[[i]]*tempfac[[i]],{i,1,iinb-iina+1}]; Print["Occupation of the initial states:"] Print[Tab,MatrixForm[tempfac]] Print["S(Kappa,Omega)=",MatrixForm[sqwtabl2]] 偏極中性子に対するS(Kappa,omega) (* 6.4.1 Definition of the Pauli matrix. *) pauli={{{0,1},{1,0}},{{0,-I},{I,0}},{{1,0},{0,-1}}}; pauli2=Table[(pauli[[i]].pauli[[j]]),{i,1,3},{j,1,3}]; (* 6.4.2 Input directions (theta,phi) of incident(1) and scattered(2) neutrons. *) th1=0; ph1=0; th2=Pi; ph2=0; nsdir1=N[{Sin[th1]Cos[ph1],Sin[th1]Sin[ph1],Cos[th1]}]; nsdir2=N[{Sin[th2]Cos[ph2],Sin[th2]Sin[ph2],Cos[th2]}]; Print["direction of incident neutron spin = ",nsdir1] Print["direction of scattered neutron spin = ",nsdir2] (* 6.4.3 Spin functions of the incident and scattered neutrons. *) sp1 = N[{Cos[th1/2]*Exp[-I*ph1], Sin[th1/2]}]; sp2 = N[{Cos[th2/2]*Exp[-I*ph2], Sin[th2/2]}]; check1=Chop[(Table[Conjugate[sp1].pauli[[i]].sp1,{i,1,3}]-nsdir1)]; check2=Chop[(Table[Conjugate[sp2].pauli[[i]].sp2,{i,1,3}]-nsdir2)]; Print["incident neutron spin state:",Chop[sp1]] Print["scattered neutron spin state:",Chop[sp2]] Print["check1=",check1,Tab,"check2=",check2] (* 6.4.4 S(q,w) for the incident neutron spin to the scattered neutron spin. *) nsig=Table[Chop[Conjugate[sp2].pauli[[i]].sp1],{i,1,3}]; Print["When polarization analysis performed:"] Print[Tab,"=",nsig] sqwsif=Table[Abs[nsig.qperp[[ii,jj]]]^2,{ii,1,Length[qperp]},{jj,1,Length[qperp[[1]]]}]; Print[Tab,"S(Kappa,Omega;Chi->Chi')=",MatrixForm[sqwsif]] (* 6.4.5 S(q,w) for the incident neutron spin state. *) nsig2=Chop[Table[Conjugate[sp1].pauli2[[i,j]].sp1,{i,1,3},{j,1,3}]]; Print["When polarization analysis not performed:"] Print[Tab,"=",MatrixForm[nsig2]] sqwsi=Table[Chop[Sum[nsig2[[i,j]]*Conjugate[qperp[[ii,jj]][[i]]]*qperp[[ii,jj]][[j]] ,{i,1,3},{j,1,3}]],{ii,1,Length[qperp]},{jj,1,Length[qperp[[1]]]}]; Print[Tab,"S(Kappa,Omega;Chi)=",MatrixForm[sqwsi]] (7) 電荷分布の視覚化 (* Spin-Orbit and Zeeman eigen states obtained in sec. 5. *) vals=sovals;vecs=sovecs2; Do[Print[i,Tab,vals[[i]],Chop[vecs[[i]]]],{i,1,Length[vecs]}] (* 7.2.1 Input number of the wavefunction you want to plot. Coefficients are calculated. *) wfnum=3; Module[{v,L,ML,S,MS,vd,Ld,MLd,Sd,MSd}, cfac=Table[Table[ Chop[Sum[ v=vecs[[wfnum,ii,3]];L=vecs[[wfnum,ii,2,1]];ML=vecs[[wfnum,ii,2,2]]; S=vecs[[wfnum,ii,2,3]];MS=vecs[[wfnum,ii,2,4]]; vd=vecs[[wfnum,jj,3]];Ld=vecs[[wfnum,jj,2,1]];MLd=vecs[[wfnum,jj,2,2]]; Sd=vecs[[wfnum,jj,2,3]];MSd=vecs[[wfnum,jj,2,4]]; If[S==Sd && MS==MSd && ML-MLd==-q, (Conjugate[vecs[[wfnum,ii,1]]]*vecs[[wfnum,jj,1]]*rmcef[k][[v,vd]] *ClebschGordan[{Ld,MLd},{k,-q},{L,ML}]/Sqrt[2*L+1]*(-1)^q *Sqrt[(2*k+1)/4/Pi]),0] ,{ii,1,Length[vecs[[wfnum]]]},{jj,1,Length[vecs[[wfnum]]]}] ], {q,k,-k,-1}],{k,0,4,2}]; ] Print["cfac=",cfac] (* 7.2.2 Function of the angular part of the probability density. *) rc[th_,ph_]:=Sum[Sum[cfac[[k/2+1,k-q+1]]*SphericalHarmonicY[k,q,th,ph] ,{q,k,-k,-1}],{k,0,2*l3d,2}] SphericalPlot3D[Evaluate[rc[theta,phi]],{theta,0,Pi,Pi/25},{phi,0,2Pi,Pi/25}, PlotRange->All, AxesLabel->{"x","y","z"},ViewPoint->{7.299, 4.552, 2.646}] ]