From 85d4664d9698533fbe35caef503b0377eb0995f5 Mon Sep 17 00:00:00 2001 From: aichhorn Date: Mon, 27 Apr 2020 16:22:11 +0200 Subject: [PATCH] started SOC tutorial --- doc/tutorials.rst | 10 + .../images_scripts/Sr2MgOsO6_SOC.indmftpr | 20 ++ .../images_scripts/Sr2MgOsO6_SOC.outdmftpr | 31 +++ .../images_scripts/Sr2MgOsO6_SOC_DOS.png | Bin 0 -> 11563 bytes doc/tutorials/sr2mgoso6_soc.rst | 214 ++++++++++++++++++ 5 files changed, 275 insertions(+) create mode 100644 doc/tutorials/images_scripts/Sr2MgOsO6_SOC.indmftpr create mode 100644 doc/tutorials/images_scripts/Sr2MgOsO6_SOC.outdmftpr create mode 100644 doc/tutorials/images_scripts/Sr2MgOsO6_SOC_DOS.png create mode 100644 doc/tutorials/sr2mgoso6_soc.rst diff --git a/doc/tutorials.rst b/doc/tutorials.rst index c24133fd..a5dd7ea8 100644 --- a/doc/tutorials.rst +++ b/doc/tutorials.rst @@ -22,6 +22,16 @@ Basis rotations: Sr2MgOsO6 without SOC tutorials/sr2mgoso6_nosoc + +Sr2MgOsO6 with SOC (non-magnetic) +--------------------------------- + +.. toctree:: + :maxdepth: 2 + + tutorials/sr2mgoso6_soc + + Full charge self consistency with Wien2k: :math:`\gamma`-Ce ----------------------------------------------------------- diff --git a/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.indmftpr b/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.indmftpr new file mode 100644 index 00000000..de341fa1 --- /dev/null +++ b/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.indmftpr @@ -0,0 +1,20 @@ +5 ! Nsort +2 1 1 4 2 ! Mult(Nsort) +3 ! lmax +complex ! choice of angular harmonics +0 0 0 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +cubic ! choice of angular harmonics +0 0 2 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +1 ! SO flag +complex ! choice of angular harmonics +0 0 0 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +complex ! choice of angular harmonics +0 0 0 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +complex +0 0 0 0 +0 0 0 0 +-0.09 0.43 diff --git a/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.outdmftpr b/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.outdmftpr new file mode 100644 index 00000000..815e1920 --- /dev/null +++ b/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.outdmftpr @@ -0,0 +1,31 @@ +Density Matrices for the Correlated States : + + Sort = 2 Atom = 3 and Orbital l = 2 +Writing the matrix as : [ block up/up | block up/dn ] with + [ block dn/up | block dn/dn ] + # For the Up/Up block : + 0.003574 0.000000 0.000000 -0.000000 -0.000000 -0.000000 -0.000000 0.000000 -0.000000 -0.000000 + 0.000000 0.000000 0.254943 -0.000000 -0.020080 0.048438 -0.000000 -0.000000 -0.000000 -0.000000 + -0.000000 0.000000 -0.020080 -0.048438 0.011530 0.000000 0.000000 -0.000000 0.000000 -0.000000 + -0.000000 -0.000000 -0.000000 0.000000 0.000000 0.000000 0.363933 -0.000000 0.157594 -0.000000 + -0.000000 0.000000 -0.000000 0.000000 0.000000 0.000000 0.157594 -0.000000 0.363933 -0.000000 + # For the Down/Down block : + 0.003574 0.000000 0.000000 0.000000 0.000000 -0.000000 -0.000000 -0.000000 0.000000 -0.000000 + 0.000000 -0.000000 0.254943 0.000000 0.020080 0.048438 -0.000000 0.000000 0.000000 -0.000000 + 0.000000 0.000000 0.020080 -0.048438 0.011530 -0.000000 -0.000000 -0.000000 0.000000 0.000000 + -0.000000 0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.363933 0.000000 -0.157594 -0.000000 + 0.000000 0.000000 0.000000 0.000000 0.000000 -0.000000 -0.157594 -0.000000 0.363933 0.000000 + # For the Up/Down block : + -0.000000 0.000000 0.000000 -0.000000 0.000000 0.000000 0.028268 0.000554 -0.028268 -0.000554 + -0.000000 0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.111512 0.018026 -0.111512 0.018026 + 0.000000 0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.014697 0.018147 0.014697 0.018147 + -0.028268 -0.000554 0.111512 -0.018026 0.014697 0.018147 -0.000000 -0.000000 0.000000 0.000000 + -0.028268 -0.000554 -0.111512 0.018026 -0.014697 -0.018147 0.000000 0.000000 -0.000000 0.000000 + # For the Down/Up block : + -0.000000 -0.000000 -0.000000 -0.000000 0.000000 -0.000000 -0.028268 0.000554 -0.028268 0.000554 + 0.000000 0.000000 -0.000000 0.000000 -0.000000 0.000000 0.111512 0.018026 -0.111512 -0.018026 + 0.000000 -0.000000 -0.000000 0.000000 -0.000000 -0.000000 0.014697 -0.018147 -0.014697 0.018147 + 0.028268 -0.000554 -0.111512 -0.018026 0.014697 -0.018147 -0.000000 0.000000 0.000000 -0.000000 + -0.028268 0.000554 -0.111512 -0.018026 0.014697 -0.018147 0.000000 -0.000000 -0.000000 -0.000000 + +The charge of the orbital is : 1.99583 diff --git a/doc/tutorials/images_scripts/Sr2MgOsO6_SOC_DOS.png b/doc/tutorials/images_scripts/Sr2MgOsO6_SOC_DOS.png new file mode 100644 index 0000000000000000000000000000000000000000..20a6f078bf681f90ce1ea6c7917149790cfa42d6 GIT binary patch literal 11563 zcmdUVXH*nz6D62Xi4vQv5(N>OC|NQ{5)qIjVaPdyL=hQ5dPH&(5D*j*L~@WULk=T3 zOO%Y{4BMc-zTckR{k7lOvz#-`%u`jjZr!Txo_VH&@2SY)U#7Z@g@uK$ATNC%3+pU^ zg>}Z{;#qhlo!9*#{3Li!>4D7g@iF}Jp9lQzOK9lu@Ngtn7uJdG_&6QQ>Dm`_9Qixn zJF!Brs;jI0{rx2+B}ufuufdU*9OSi~v9NGSP9E4;@vq@DSTZg$4_zKuy4aXvaay`K za4Or|JDbX?a~fG9EG+Du>^KodHV9W6BNt11J7-QUOH(^;Yfg11Ze?>-XH|YqX;pOt zbycZbM&?cePZd>QK*AG130D_580!znbXfQW*hF?>QnxpAaW`@@#gcO}LYQ*!-2!iM z^S$)cIjIZp#6jBB8R2BeRy!^g+BKl=LKQC~LPFJ+RgG?%wP+z2I4;9j&Nyzc2{{I+CbAyVeK1|f@ zt%U=r+*W^BYKt9ZD# z?{6q-J8MpszchbX9hPd)9G_1$<|*h~K?n6;|kz;pAvJccuH;iGS5Ti*&j%Ev6pAZBcgnd{W=;tR8# zcd71m52}(wt&a^{HxKu-Xn!7NCB^oye8Ts2Toi(7Jxw)&B&ZgC42nOVKpAnY?eA_y z@9Z6Xa1_NuG-4s}zG;dN^&wXrY;EwCAOm43)FJD;+dPF!cfmiwN{15$rk5|W9DDn2 zzQb9koLE<*r|)LL*<|ELZ~xfVt+?O6P<2dcZh=;{g-bjqQwId_Qc6{^+_!(_+nSJO ze@K-&#>IrbSGgL%)|SlEibJ?L?7L+C$G4zI)i0Mnna;QE-G$9mHFKy+b%LHgfASDA zfB58$1%P-}w3mr4N=_Bun&_a4LU(dfnZ32udI9T9M2de>K~)PN-V2OXnel?iy?HUj zd+YT{U9U&B<@f&$krwy#37vR8(Y74lD^-yT#Y zfEyu5TCEE;>}ZD50gRt_=aKzsAkbOIg=TO{B|PkTg0jdm*kO6RC(ydY1*rwThBHSN z?FIUUzdMe&Sfgrh1j$WF*8yp8Ejtk+uBIrXBPk*LOsWhspd%Mc5W|IL&E70OVq96b z=Iox8J1pTpEJo2}muH;x{1okv#-}iUG``B>nF97MwKm6>1E(GMWi1UkJ>2f7YzQGf z4#marxr3+2SYz;`zN7`Wnfe3l!rPxAvl9_ic;$lk4F&Z(+rsK&NGfuKji5prZEr!$ zLKzAAVaxxG#gD>KoSIAcm5g4{6+ui(y!ukPmHkhq0R+?-j08uU!wO?a@DX347_2#2 zsL&qla=njKc%fCzmme$O(kHWNw7pIg(iP0xa)M6)#ZD$m)3BKK(~W~I|2=>WHpHwC z#S%X02_0C-5}aJ*Uxr5{#5me*8DgI;iBju zI8TLTk$!>TPWU1y`&F8ApFrGy75~GEZL)$uyl$O&7&Coqr? ze;fqI7MYp(FlQt^hx7bheq@4Iyb#_-X#CqY+X?WfUaPfbDwz&z{n+ufd!2}+Vi<7p zqHFA4Q&JKo1ToEkU4Gmf32do=(Wr}hDerlI0r z9ysV&&&%#PiBnr5VVQE=;HJn|X;mIskP=dp+u@ONn=wraI$>r02SE=VdynM^I&(tX z?E*UJ3Q5H%P#?^RbMjKO1(-IG$Q4%;i>Ah^tp2gHGU!H8Y=J*m zr%pyt;wrj?rWVIb@)DezS5}M~Qo96QVSqL42Q;LCN#qFuL{8m}=zNl`XT&F&`I)(8 zGvCtYGDk=%Z4jUvRlNndWS_>)lA(hncZ2P5(DJZ0XdwMvwR0RcI`_RKSD(M^D4@V7 zaOHh-D;M&i`|2rvtTvkK6oeQ047wv21F4_!#gU%0SM+KK7J@M-xDu+--S#eG3{Ti$ zo`?tTBktvUj<+(Ma>9idj1fPhU+}~LCovpB*2Ygj;;+&n)OD|a`b=$YM1AE0tAR_# zsVnI3L;MLhMMqn|YDvW4p75ik3dmpHQ1DaZmzR4aHV1_?KZ&o32wWd@vo(_GiVR#yRdraX8(d3XfY?DE{S_ zm*Jq&RT{&mutfJ$S*REyIB1pS9^pW0A(3!babKW2!-lX3UzdJ;Ma*Fb?JREI6Il~7 zgmm6Yb8^XwNa9zZS^y)CVIE*o@-$65u&tx#_5`H44+~D)i6TUFUuA+-!`8C-rvF;k>N2Vtt$mf(ld9t9Z$VbwFU=S_=EFh@sVui<7tK|-)MfU>|;C=Dj({AU>sH4As6o%AV z8vMlbHp>E?C}`krA&TIdD8$IntPo!HX;Uy9 zP#Zq=y=oe1Utug=&652a(!r0c>~2~%SD0}Vt|Cqy_GQFZlN#a52uI}2?xMK}we}1} z(8|6{Uw!sd?nBt<2Tt$XH^0_68bbP4mxl& zpAw>5+YeFTutrbWfPC)a8V|DD^ zmvN~9ZYa8LgnK9x+97jd94Hs3%$mH1U>$h&Kl?B|*5tc~cJLo^fJTy`9tnZ|EfduE zJo$NUM8!!DGG^acr{QeMCu!hoBlO-Wyery|_$)MgVJHa#_2S+Ge82e+l<|As3*FGx zz(N#+N7%|`Uj0u?T%9K~pX?f>V!O?PaGeseadygEVkLWv{#|tYS^BYv8biP5(a^X` zywYEhbY@-^_q+c9}Z6B9fGZzzwysE1bcwhm(<~`N*8GO@1UD%#CEtZajU!B zb03Nrf9U?a7CxYap=7U|I{9lj``=p_8 zaSWIx%wD(^x-Pw?%gJxtS|Q<*jhiBv{E@W4b-(d9bU{|OHOr<$qD+xrGQfXw+YwytjQ2ZH}x(%J|=$^;ZgpckxtMm+t z_+6@nmTsTQe9g^?o(`%l4Ij{6VZkovui^`Rg}ov&w>R=N6*21PkhpLn{o z5)1vfWW};CGun1NRn8S?q9rHYM=1Q2st7*y%P~B-YMDT8^Ryy(h0B1@ooDKtQ!Zhp z($iMbC2L&7MdC3dgQ3gh*qI_TyXLu=fE{oVLG0i7dTYosM44 zVs4#`*J)Pn*=90yuHsPdv#B8q`s*k{r*-%4$3mnw!HjXNLZV4;*l|$fnnmec`|iv8 z8dhmOLA+v(9K&j6C1;DI!<9%A=|U})yaWv)c30iY<$gn{Pg3e70o0?Z06jh)9wwe) zY~AY$Fh@3P&H(i%Aaly5{`KUk3v$F^l*>;z!Wc?%fJgby{zD^`Cy#K*ayrGQOj2e~8mQ=zUZRtYdzq5kK36=UyExA0C_<&^*}U*943xFsHfdPtcjx0=&OS z9=`&emI`RH7l^l6NziX#dEPJfC~L`m;+5_C5MsEnOp}i(vopPQ7SRdoqUlD0;LRz! zSkzwF<+LoY^O5*8B>N@=Chy~mWjzBem|v){@2>z1Lnr9=6)J@wVS#>4ixk|uU{?Y;yCkMJb~r*=L4NtISHW7mN;rK>9JD4^e_ zrXdo{<$Zc8>bodYE~J9d5&TiM$d#aYGs;qY5g*yH2Rx*zIjyp*8JRddppyc|QxC5U zl(df*Eo*WDX+9|Ck}YNO;LAEySF*_l9F^#T>Q6AQ(YBgevg0cy|H#0PgV5If?sO#8gwUNxPr-2Rpv6;r=P86{yk18QaqZ_>q3YeL!Dd z4)J-sTUCFOOr~AW=zJeod2Q*=(%=&T%Cr?eJ`!UICG1 zK#xG#Xx$$UBT7mst3I=uH~Eq6T@l9oHHWSmtp5)0co8d^p9~&!Rg~~<{>!B98VE8w z^T-0ThPV9#PLrv-PV=H|C25bpfP2m1{u~siGG@S%b9L!(SH8frYe@N#_4>OIQjsGw+ z!`QGy1s1;{2hz^}buSg%`{Hv?Gys%0ge;Xv^woQYX>0?}LMb?dB&&`B=nUkEdBnEo3d zj2h>3!nN`Tb-p(@M9M~6rdof9Pe^MRiQ895=0II}_%LaIW?1Oo$me`0l=eRI$y#*Z zxf<}nF7mCRvmLG8?jJ#%ou5|XGJ^%*W5JwH{&y~DjOg7zP*?ed5Iy^0n<7roEgQ9` zfX8wBn~Yse4{9$=LPH>i^f;0v$NNMzrm_$3;}qO>jsnxq+v|77wzO51sflyHpPvbz zNOr1?3cZ24nJ=wH#>dk{FjkdYU|3Sg7nAVGbk_Q?@y!S7LdRUeVuD{#*Mx*ViTuH; z+srN4fB*KE`&X%lo%tU0fM)|74<@Uyi2F zci(}?a-N;kc2#K4{N0&2ZKJXC%36spbj> z<@+_}YFBF3v%!>tr_>!%6a?8uOWLOl+%a-?lq`n@uZA_e-8_^XA3~8Ne3(x4AzObLK6+yj@T`}HzNUazdz!c1Tno@)J zPfGW1(QO^Q+Z1{#Hmx|zn?^Y=kt|G##FjfU33i#y3322b`X2m++SZmWU6q~N6vqYDlo;R<{|NfDRlL&bo04RUp~GuJmn1gPGnn!9&JE?d3oSPGEA5*5NPA2C5Zc3Mai0BGYgnY_-FLo>)f* zfa60g?hAq2$5F(?!OE^S~gyW*)e0KOnfsbb%2K9zJYom|vQWW>8!pT71Hw{p6Ss~op9WeYhTUy`}ah3(T) zy0PndyBVVl{IM&$ZA_s$3&`p!&|(&^k!D7EK9UN6Rm)TBx!FHLPM zp~AehTk%w{v4)q%S3b5oi-45wslTwr&7U>)zu!2A12DX#a+>5R+0wDWQ%5Ov-h*wQ zW^FoHkM=p~QWi=}!APEo_tyl(0#J|<-Q4>OqCtvBnDn&+YV^)v9pI1TfWktoE~og6 zlj*Y&LN)Ta8Zzd2uLj;2?@L-3c*w*wuM+~3RqoUQ!zrgr4`s})%y%iEl%tr6mf@Xq zl55q(-$hNlY?$hYuonD9FcFlSz*kGlI(6dj5 zzzuq{uMt#hsN4upQh-rH&d%GqMA|%SB6I1jGAe~m3Ng!aWph`yzXjOXM8N=D!R=Kq zJq_`g4e>6;N}Go`T>7(2iH<2HSbYp=yB8BxK^!SY)y-bSxHvk`M{S zVY!k(jP|&$ORDn$_Cgh=cE3_;GqKfcBfTb-j6n`UzbEOzG3%_X6mp@YotTS0_=B59#yewOV2>3EMGTl{ETd z5>W?7+<4nAFfq<-6`i##I0WPzSAefY4@P(y7oRmBKRytDZo)WxV}Os+a4tv2AtF!j{$#-1EFqZ`Kv*&iIel^@OExrn;*@(0ookLtLI zmmf)~95dP+`&y8vr-=b^?l^UvuGReK{4<|dqU01!G2WLSokiLw^Zg4T1onRKG#sD8uFrzX6MhXDtO!;<#H{0TIH=Q^*qw zFE!W+*tk@W=6biK3e@S{;zJIK@ylBDUV?uSwpNr?eIntaL(tN8F$VMwdgv)Wr4OY% zz@RTy309#2&+-Sjd>L}#7!pIG#OOL=kE@=n0E1_At6@5K`PO+Oyd^vV1<}kS-m&vx zBk;Er?g5}^K%GNcuD5U{LZ0eRsvD>o~}ADAf#@QwlWcUvex ziwhGN1{gcC&ZtAKYMd*n)_T!V~ERRHF=jarQ>>AMFJQ; zSzu6CqYh9cCdCADBkeyx@L@3oNbMG)0Rt{sR*YYQqDyjh$?V3Ln{#vLL|vo7+Ju2Y zKvP$9+~hi zBZjN)o<-<=ZJ54B@MFtI%g-PD5Uj3CnF4&endQLkkN!S2(_?bIE_HI<>D7`tzwG=O zVDkE1@2-yt(o%V{H@&%^Q!>!KjkADsZ#<=#>csC(zsZGU z=+f!Hi2wq)vH9SrBYpS_Mij?I!e`g5&pS^(7`-f<&CEkdmN@tEM*xW9$%7{Qf)0{U zv07U(gJ#!p3)o^Lk_ocN4) z1|8I3iTLQ-^aT92Stw^rFu>IKg}*b+yL-x>TXw0{>?1f~Bu;KCs-^&6Snlv@ZU#%E zOX_LjjQFkA1c0=onFO`FO?KPppy$D*>;53`-1cr+v06L=Qu9+(@XAYJS2pVa7-)%V zNH<0)$oNJ>K}_Faf+~5ugXAUA%^Dy-bD(JaJmQgnnE}^gwB=3OuxgkwxtL%eh?6V1PR-zg94cM;Naezw&4tXlpA5w}l z6z)=qi$A@F`H-!nxQIdLwg(|X3qD=Eq$$1vU(b*5*(*x&T2?~R_(BC|v)L(izP!kr zb$URno%v?EXRE+@te{CkRdA4qnZGKO(LMIoMUTCeACYM${kzxA@y4d&a+zLQaq7*w z>tG_Ecx)y*ZQS_6Vv3ob)f%86%xZJY!gZn(oV_7OUHf%rJ8b8C+UxyKH4`b+_4OwG zm23<4*Bo2Fw_hbGV#HLF4Z74sM4D`({$RX6vb!6Fx&J6THL#~61h@Nn9BNgd0^q!m znfq9c-&7?{pn2S>J4v>eX1Zs3R^3CjIh8(Dp=6wBOQ7Po27}_*2Zz%6p2M zc{v1QR@eD!t?1mMT-D{zu}&p4__aAU0?k>!?}ZbWri8Ar-UYbi%VEQ zMCA%s7n^`NjtpGMZ@U7`Azcpdg;1k6p0{9bdw%rEs_kq@-ONC-_*~nNs8djd-Y|0*}m!p`wt*J>z zc?VmZbdM``x<0CWu(ns9TE3u*EFJjy8P!p(o9S=w_aeEyeapoMCm%;N*V(?bvKt3l z4o4m8AonixZ`-y$_)t}&30grb{@blp$;ruMqs0sJLo4!LzH#ltFk0lzW1vgjKH>7KeuhaMU3*(sl-o})0Pa;Oi^_MbPsMwff`*7cHq zc0sEMzk^M(4dSd6wl@OTW#8X?A#v7->g8D@V3ermmy?unAhdluxIBLx}$3^?9kxTmf0z;KEEIP(^Sjd2RGlm zd&n?k=Vx2cqmC?2YA;%FFh&+0}GV0D$RZ-QKay(xCV`)t0ctU(&q;l5a zFwJ*khp*G{U(<*6YFzgasr6$8AN%JgzDu*7U>``A2V)#v_6cIi>Z-U~tC?fsKv3Nk9vIg*A?{|At^mMZ`N literal 0 HcmV?d00001 diff --git a/doc/tutorials/sr2mgoso6_soc.rst b/doc/tutorials/sr2mgoso6_soc.rst new file mode 100644 index 00000000..596e7d29 --- /dev/null +++ b/doc/tutorials/sr2mgoso6_soc.rst @@ -0,0 +1,214 @@ +.. _Sr2MgOsO6_SOC: + +In this tutorial we will discuss how to set up and perform a calculation including spin-orbit coupling. As an example, we again take Sr2MgOsO6, as we used for the discussion of basis rotations (:ref:`Sr2MgOsO6_noSOC`). + + +The full script for this calculation is also provided here (:download:`Sr2MgOsO6_SOC.py `). + +DFT (Wien2k) and Wannier orbitals +================================= + +DFT setup +--------- + +First, we do the DFT calculation, using the Wien2k package. As main input file we have to provide the struct file :file:`Sr2MgOs6_noSOC.struct`: + +.. literalinclude:: images_scripts/Sr2MgOsO6_noSOC.struct + +The struct file is the same as for non SOC calculations. However, the DFT calculation is done a bit differently now. First, we need to do a spin-polarised DFT calculation, but with magnetic moment set to zero. In the Wien2k initialisation procedure, one can choose for the option -nom when ``lstart`` is called. This means that the charge densities are initialised without magnetic splitting. After initialisation, do a full SC DFT run:: + + runsp + +After this SC cycled finished, we continue with the spin-orbit calculation. We initialise it with the Wien2k command:: + + init_so + +In this initialisation procedure, be sure to select a spin-polarised calculation. Then the SC cycle including SOC is done by:: + + runsp -so + +After the SC cycled finished, you can calculate the DOS with SOC iuncluded. It should look like what you can see in this figure: + +.. image:: images_scripts/Sr2MgOsO6_SOC_DOS.png + :width: 700 + :align: center + +Wannier orbitals +---------------- + +As a next step, we calculate localised orbitals for the d orbitals. We use this input file for :program:`dmftproj`: + +.. literalinclude:: images_scripts/Sr2MgOsO6_SOC.indmftpr + +Note that, due to the distortions in the crystal structure and the SOC, we include all five d orbitals in the calculation (line 8 in the input file above). The projection window is set such that all d orbitals are included. + +To prepare the input data for :program:`dmftproj` we execute lapw2 with the following options:: + + x lapw2 -up -c -so -almd + x lapw2 -dn -c -so -almd + +Then :program:`dmftproj` is executed for SOC:: + + dmftproj -sp -so + +At the end of the run you see the density matrix in Wannier space: + +.. literalinclude:: images_scripts/Sr2MgOsO6_SOC.outdmftpr + +As you can see, there are a lot of off-diagonal elements now, in particular also off-diagonal in spin space. This is just telling us that spin is not a good quantum number any more in the presence of SOC. + +We convert the output to the hdf5 archive, using +the python module :class:`Wien2kConverter `. A simple python script doing this is:: + + from triqs_dft_tools.converters.wien2k_converter import * + Converter = Wien2kConverter(filename = "Sr2MgOsO6_SOC") + Converter.convert_dft_input() + +This reads all the data, and stores everything that is necessary for the DMFT calculation in the file :file:`Sr2MgOsO6_SOC.h5`. + +The DMFT calculation +==================== + +[TBC] + +Rotating the basis +------------------ + +Before starting the DMFT calculation it is beneficial to look a bit more closely at the block structure of the problem. Eventually, we want to use a basis that is as diagonal as possible, and include only the partially filled orbitals in the correlated problem. All this can be done using the functionalities of the :class:`BlockStructure ` class, see section :ref:`blockstructure`. + +We first initialise the SumK class:: + + from triqs_dft_tools.sumk_dft import * + SK = SumkDFT(hdf_file='Sr2MgOsO6_SOC.h5',use_dft_blocks=True) + +The flag *use_dft_blocks=True* determines, as usual, the smallest possible blocks with non-zero entries, and initialises them as *solver* block structure. In order to disentangle the :math:`d_{x^2-y^2}` and the :math:`d_{xy}` orbitals, and pick out only the partially filled one, we do a transformation into a basis where the local Hamiltonian is diagonal:: + + mat = SK.calculate_diagonalization_matrix(prop_to_be_diagonal='eal',calc_in_solver_blocks=True) + +We can look at the diagonalisation matrix, it is:: + + >>> print mat[0]['down'] + [[ 1. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j ] + [ 0. +0.j -0.985+0.j 0. -0.173j 0. +0.j 0. +0.j ] + [ 0. +0.j 0.173+0.j 0. -0.985j 0. +0.j 0. +0.j ] + [ 0. +0.j 0. +0.j 0. +0.j 1. +0.j 0. +0.j ] + [ 0. +0.j 0. +0.j 0. +0.j 0. +0.j 1. +0.j ]] + >>> + +This transformation is already stored in the SK.block_structure class. The next step is actually not needed for a DMFT calculation, but it is good to do this check to see what the transformation does to the local Hamiltonian. We can calculate it before rotation, rotate it, and look at the 2x2 block with and without off-diagonals:: + + eal = SK.eff_atomic_levels() + eal2 = SK.block_structure.convert_matrix(eal[0],space_from='sumk', space_to='solver') + + print eal[0]['up'][1:3,1:3] # prints the 2x2 block with offiagonals + [[ 0.391+0.j -0. -0.815j] + [-0. +0.815j 4.884-0.j ]] + + print eal2['up_1'] # prints the 2x2 block after rotation + [[0.247-0.j 0. +0.j] + [0. -0.j 5.028+0.j]] + +So the local Hamiltonian has been diagonalised. From the other entries we can see that the *up_0* block and the [1,1] entry of the *up_1* block correspond to :math:`e_g`-like orbitals, and the others are the +:math:`t_{2g}` orbitals that we want to keep. So we pick the according orbitals in the block structure:: + + SK.block_structure.pick_gf_struct_solver([{'up_1': [0],'up_2': [0],'up_3': [0],'down_1': [0],'down_2': [0],'down_3': [0]}]) + +We can now look at the final result:: + + print SK.block_structure.convert_matrix(eal[0],space_from='sumk',space_to='solver') + {'up_2': array([[0.156-0.j]]), 'up_3': array([[0.156-0.j]]), 'up_1': array([[0.247-0.j]]), 'down_3': array([[0.156-0.j]]), 'down_2': array([[0.156-0.j]]), 'down_1': array([[0.247-0.j]])} + +We see that we arrived at a structure with 3 orbitals per spin only, and blocks of size 1x1. + +The interaction Hamiltonian +--------------------------- + +We now set up the interaction Hamiltonian. Since we want to rotate the interaction matrix into the local basis, we are using the Slater convention for it:: + + from pytriqs.operators.util import * + from pytriqs.operators.util.U_matrix import * + + U = 2.0 + J = 0.2 + U_mat = U_matrix(l=2,U_int=U,J_hund=J,basis='other', T=SK.T[0].conjugate()) + +In the last line we use the Wien2k convention to write the U matrix in the cubic harmonics. Next, we want to rotate it, and pick out only the relevant :math:`t_{2g}` orbitals. At this step we need to give the indices of the wanted orbitals in the *sumk* basis:: + + U_trans = transform_U_matrix(U_mat, SK.block_structure.transformation[0]['up'].conjugate()) + indices_to_pick_sumk = [1,3,4] # these are the relevant indices in the sumk basis + U_red = subarray(U_trans,len(U_trans.shape)*[indices_to_pick_sumk]) + h_int = h_int_slater(['up','down'], indices_to_pick_sumk, U_red, map_operator_structure=SK.block_structure.sumk_to_solver[0], off_diag=False, complex=False) + +Now we have the interaction Hamiltonian for the solver, which we set up next:: + + from triqs_cthyb import * + import pytriqs.utility.mpi as mpi + + beta = 40.0 + S = Solver(beta=beta, gf_struct=SK.block_structure.gf_struct_solver_list[0]) + + # Solver parameters: + p = {} + # solver + p["random_seed"] = 123 * mpi.rank + 567 + p["length_cycle"] = 200 + p["n_warmup_cycles"] = 100000 + p["n_cycles"] = 1000000 + # tail fit + p["perform_tail_fit"] = True + p["fit_max_moment"] = 4 + p["fit_min_n"] = 40 + p["fit_max_n"] = 100 + + +The DMFT loop with automatic basis rotations +-------------------------------------------- + +After these initialisation steps, the formal DMFT cycle is very similar to a calculation without these basis rotations, since these rotations are done automatically, once the :class:`BlockStructure` property *transformation* is set, see :ref:`basisrotation`. + +The DMFT loop itself looks very much the same as in :ref:`SrVO3`:: + + # double counting correction: + dc_type = 0 # FLL + # DMFT loops: + n_loops = 1 + + for iteration_number in range(1,n_loops+1): + + mpi.report("Iteration = %s"%iteration_number) + + SK.symm_deg_gf(S.Sigma_iw) # symmetrizing Sigma + SK.set_Sigma([ S.Sigma_iw ]) # put Sigma into the SumK class + chemical_potential = SK.calc_mu( precision = 0.01 ) # find the chemical potential for given density + S.G_iw << SK.extract_G_loc()[0] + + if (iteration_number==1): + # Put Hartree energy on Re Sigma + dm = S.G_iw.density() + SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) + S.Sigma_iw << SK.block_structure.convert_matrix(SK.dc_imp[0],space_from='sumk',space_to='solver')['up_1'][0,0] + + # Calculate new G0_iw to input into the solver: + S.G0_iw << S.Sigma_iw + inverse(S.G_iw) + S.G0_iw << inverse(S.G0_iw) + + # Solve the impurity problem: + S.solve(h_int=h_int, **p) + + # Solved. Now do post-solution stuff: + # Set the new double counting: + dm = S.G_iw.density() # compute the density matrix of the impurity problem + SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) + + # Save stuff into the user_data group of hdf5 archive in case of rerun: + SK.save(['chemical_potential','dc_imp','dc_energ']) + +The only difference to the other example is in the initialisation of the real part of the self energy. We cannot just take an element of the *dc_imp* array, since this array is stored in the *sumk* structure. Therefore, we first need to transform this matrix into *solver* space, and then take the appropriate matrix element. After the first iteration (here done with 24e6 MC sweeps), you should get self energies like this: + +.. image:: images_scripts/Sr2MgOsO6_noSOC_Sigmas.png + :width: 600 + :align: center + +The two :math:`d_{xz}` and :math:`d_{yz}` orbitals are degenerate (blocks *up_2* and *up_3*), whereas the :math:`d_{xy}`-like orbital is different. + +A complete python script for this tutorial, including some more input/output, is available (:download:`Sr2MgOsO6_noSOC.py `). When running the script, you will encounter warnings during the transformation from the *sumk* to the *solver* basis. These warnings just reflect that the off-diagonal elements of the full Greens function are not zero at all frequencies, although the local Hamiltonian is. In that sense, we still do an approximation when restricting ourselves to the :math:`t_{2g}`-like orbitals. \ No newline at end of file