|
平均場近似 （self-consistentな固有状態の計算）
基本例：
- – Note –
- モデル：AB部分格子間に反強磁性相互作用が働き，z軸方向の磁気モーメントが互いに逆向きに発生して秩序化している．x成分とy成分はゼロとする．さらに，z軸方向に磁場がかかっており（平行磁場），AB両サイトの磁気モーメントが異なる絶対値を持つ（磁場方向のモーメントは伸び，反対方向のモーメントは縮む）．
calc: 初期値を入れると計算結果を返すモジュール．ここでモデルを定義．
jzai: A格子の<Jz>の初期値，jzbi: B格子の<Jz>の初期値．
jzaf: A格子の<Jz>の計算結果，jzbf: B格子の<Jz>の計算結果．
jab: 交換相互作用定数．
eps: 収束判断しきい値． |(計算結果ー初期値)/初期値| < eps なら収束とみなす．
-
(* Module for the mean-field calculation, a two-sublattice model with AF interaction between Jz *)
calc[jzai_,jzbi_]:=Module[{jzaf,jzbf},
da=(dcef+g*0.67171*hhh*eljz)-jab*(jzbi*eljz);
db=(dcef+g*0.67171*hhh*eljz)-jab*(jzai*eljz);
{valsa,vecsa}=Chop[Eigensystem[N[da]]];
{valsb,vecsb}=Chop[Eigensystem[N[db]]];
za=Sum[Exp[-valsa[[i]]/temp],{i,1,2*J+1}];
zb=Sum[Exp[-valsb[[i]]/temp],{i,1,2*J+1}];
{jzaf=Chop[Sum[(Conjugate[vecsa[[i]]].eljz.vecsa[[i]])*Exp[-valsa[[i]]/temp]/za,{i,1,2*J+1}]],
jzbf=Chop[Sum[(Conjugate[vecsb[[i]]].eljz.vecsb[[i]])*Exp[-valsb[[i]]/temp]/zb,{i,1,2*J+1}]]}
]
err[]:=Module[{},{jzaf,jzbf}=calc[jzai,jzbi];
error=Chop[Abs[({jzaf,jzbf}-{jzai,jzbi})/{jzai,jzbi}],diff]]
(* mean-field parameters *)
eps=0.0001; diff=eps/10; jab=-5;
hhh=0.01;
temp=0.2;
(* initial values *)
jzai=-0.5;
jzbi=0.5;
(* one cycle of the mean-field calculation *)
{jzaf,jzbf}=calc[jzai,jzbi]
err[]
(* a simple method to get to a self-consistent solution *)
findans[]:=Module[{},ia=0;err[];
While[(Max[error]>eps),ia++;
{jzai,jzbi}=Chop[{jzai,jzbi}+(calc[jzai,jzbi]-{jzai,jzbi})*0.8];
err[]];
]
findans[];
Print[{jzai, jzbi}, Tab, ia, Tab, error]
例：温度変化の計算
-
(* calculate the temperature dependence *)
ans={};
temp=0.2; step=0.2;
While[temp<30,findans[];
Print[temp,Tab,Chop[{jzai,jzbi}],Tab,ia,Tab,error];
z = za*zb*Exp[-jab*(jzai*jzbi)/temp];
fe = -temp*Log[z]/2;
AppendTo[ans,Chop[{temp,jzai,jzbi,fe}]];
temp+=step];
(* Plot the antiferromagnetic and ferromagnetic component *)
ansexp=Sort[ans];
jzAF=Table[{ansexp[[i,1]],-g*(ansexp[[i,2]]-ansexp[[i,3]])/2},{i,1,Length[ansexp]}];
jzF=Table[{ansexp[[i,1]],-g*(ansexp[[i,2]]+ansexp[[i,3]])/2},{i,1,Length[ansexp]}];
ListPlot[jzAF,PlotRange->All]
ListPlot[jzF,PlotRange->All]
 (Fig. 5)
（左）AF成分，（中）Ferro成分，（右）Ferro成分を帯磁率にしたもの．点線は相互作用なし の場合(Fig. 2)．
Ce3+: Γ7(0)-Γ8(42 K), H=0.01 T.
-
(*********** advanced calculation ************) (* rapid and safe estimate of the next solution, set frac *)
ans={};
temp=0.2; step=0.2; frac=1/2;
While[temp<30,findans[];
Print[temp,Tab,Chop[{jzai,jzbi}],ia,error];
z = za*zb*Exp[-jab*(jzai*jzbi)/temp];
fe = -temp*Log[z]/2;
AppendTo[ans,Chop[{temp,jzai,jzbi,fe}]];
nn = Length[ans];
If[nn >= 2, {dum, jzai, jzbi, dum} = ans[[nn]] + (ans[[nn]] - ans[[nn - 1]])*frac];
temp+=step];
-
(* calculation of specific heat and entropy *)
ndata=Length[ans];
fengy=Table[{ans[[ii,1]],ans[[ii,4]]},{ii,1,ndata}];
ery=Table[{(fengy[[i,1]]+fengy[[i+1,1]])/2,
-8.31441*(fengy[[i+1,2]]-fengy[[i,2]])/(fengy[[i+1,1]]-fengy[[i,1]])},{i,1,ndata-1}];
ndata=Length[ery];
spch=Table[{t=(ery[[i,1]]+ery[[i+1,1]])/2,
t*(ery[[i+1,2]]-ery[[i,2]])/(ery[[i+1,1]]-ery[[i,1]])},{i,1,ndata-1}];
 (Fig. 6)
|
