From 283486eaec8bc459ee8339a94fd29ce98af993bf Mon Sep 17 00:00:00 2001 From: Sukhi Singh Date: Tue, 12 Jan 2021 16:29:42 +1100 Subject: [PATCH] SU(2) Exact Diagonalization Added a function to compute the low energy eigenstates of a local, SU(2)-symmetric Hamiltonian, in a way that exploits the SU(2) symmetry. Also added some generic SU(2) functionality e.g. functions to compute Clebsch-Gordan coefficients, Fusion Rules, Fusing two indices etc. --- .../FT_Su2CoarseGrain2To1Site.npy | Bin 0 -> 19569 bytes examples/Su2ExactDiagonalization/Su2Engine.py | 263 +++++++++++++++++ .../Su2ExactDiagonalization.py | 277 ++++++++++++++++++ examples/Su2ExactDiagonalization/Su2Tensor.py | 90 ++++++ examples/Su2ExactDiagonalization/main.py | 113 +++++++ .../spinOneFactorTable.npy | Bin 0 -> 19569 bytes 6 files changed, 743 insertions(+) create mode 100644 examples/Su2ExactDiagonalization/FT_Su2CoarseGrain2To1Site.npy create mode 100644 examples/Su2ExactDiagonalization/Su2Engine.py create mode 100644 examples/Su2ExactDiagonalization/Su2ExactDiagonalization.py create mode 100644 examples/Su2ExactDiagonalization/Su2Tensor.py create mode 100644 examples/Su2ExactDiagonalization/main.py create mode 100644 examples/Su2ExactDiagonalization/spinOneFactorTable.npy diff --git a/examples/Su2ExactDiagonalization/FT_Su2CoarseGrain2To1Site.npy b/examples/Su2ExactDiagonalization/FT_Su2CoarseGrain2To1Site.npy new file mode 100644 index 0000000000000000000000000000000000000000..cc16d159de51358c1b0ba4ed74bef32e3ca3008d GIT binary patch literal 19569 zcmbuGcl_2<|Hs=UJJ}(#GP1JAWky9fLflr1GRk>(h3L4=6K;Etgl;o?WyP1wjF1^B zdy|nU>-T(puIrqi%j5I=>-W9gZ|~vt?3~Z_xjyH*>ea=4_t|~V{hK#Cw%LhY1`ioI z$hvIVxyzV6x~$u|%V8s|Y6FHJGQtKAQM>1Wqec#)KXT}RV}?+7?Y3^`(E0!O)Bfa^ zgNBbf`k2w(2aT{H-H#r1RCU+@vjL;qA7VoWjTk<%YNG~K+xup<|Em4)!L_Mx-ly3a zzQvBs2l|#{eJkf%ckR=>Tm3&8TpfMPkoLY!{}#=fHQQrz-}bCC$NF~8FVwY9%l^&j zjq&aK_U+sIFMf*9r&9F`4?OE2zewN|`a_9~LYBH~63!{EYE#Ac zioUn=eG1`j!S7xZ)~@E??xwy;m|w>|qTkc`y$YcVe(y3NOYN70F1;N4M8B`|`xV0e z!S7!tWT^)v;p5*AT>RX#-)r*!ivB?74=RKQ2Y*PTFy8=&CSiUZ2Sh*6`9XznaPULQ zge>*2B>Y`(%D=ef4v&7Q^TP_^5y2l>CS<8cCEB;ifF_5b38N!$PzMt_m>7Z<`yg1@v;SpRja{<0+e zv-lN1ulUQOzry(|3*lA4UtK2TR~(mwx7L2eioYiM@y=gc2(JtNdI=jYjNZysKOqT! zOBd(A(cj?wjfL>0;BQXC`s2BH0wBC434cw(Tce-o{B4Es_TcYG!rJcg9`8)T+8%j1 z?u!0y=kFTe=!Mb zS2h2NFC}5VI4?&(&G}af;j6*FmV`A){S~eH=}CC`mHL@l@vle!hVyS0!ncBdyG+PZ z-$}xE=e;qZUwpO6&xrnA=ie)Y?*~7#OxUvOXC>izo$p>~$hogq{OssIaQ?$W_)+j5 zHwyFH>XRg#`jUR$ulP@+|IGQ%3*i^Re_1AEsb3}G&Q_Uw(o|05ucQCQ`ELv1cfrpo z6SCCrlkoKo&OC3~elwf=57GbV{7;2&Zty=h3iHb`FA0Ynrk}Se{+H-~b$)&!ToC+k zWkQzvdlELgdv2$WU$e(QqW{zRza+F~e2NOzoLDArS+y375}|f^pQ>0(n$cPjowb%$ zs5T_Q+Nvx*#?+UdKTKIWMhTTZkk~>rqqQeGTUc7579k0?sLB#l^QtYzD528PW{cB| zwgl1HlF|y*fh5>cDoao;s@9QFLZy!owlvLX%MhI{E3Ht=kpx>_Wdo{qJ2k7?3XBpe zeT=acX-4ZrbheVTLaj^^Y!#IyC^p!cQ9`8y(N?7yZ8f5^)uk0`4U%AMsw_dV!L=AA zR66Ob3(aUMhR7X2K3`9wgt^-TN0gZC9P0flLXsF zWgV4n(&9C&+O~`mDxH{Cp&4yEqO0~PU^}WTL9wcx7$sCXLv3f8(Rvb{?INvE zyOIR!rLqB)-uu>@Q9`9t)%wtkwj0se?$Qd?mn7I8DoaqTYEMQ9Mas{s72As}nj<>f zn^>Uwkp$aEWgXSBIB;vL0<|xrgsPoNb1JqU&1m{==8P|+HPitl!TzNZREvL}s8u_V zQ9`AY(GH>+?O>v_L!{LP4+Xr~gLohGeNr;`LbLuCUh9m*CNB~&_IZ7j`b zXA+&AC9P0rlLY&Z$`TYCJcm(2?Okl}fr_0=GunAXXXi^R)CDBLE>u~9VuKekN~q_H z4Nk4t#WbT`LUeYiv_f4*66|u74XAuyUBM`!(xGfu(u{T$(b?6~3N?--*flCkP;78K zqlCJn_!jd~#jd3p?K+~f>!lTH0!gs{sw_dV!5bJQRP7aHgE!KQb`#Or&C&{W3rVnB zRn}3pCqwNOWK|OxCDgv@zPgQOwA+c!?vPffJ4u4wrLqLYs_tf#P=_!0b7oj)R+HUB zGupjGXZJ}f)cqvE9#GkU%FnAwj1sDUv8qoi_8`q@4-uU`EUi$FkOX^FWeJK^O=gr( z!)mLl*kd%KJx+AiB&|?SkOX^DC8(D9(U(^>g;7G)o($ht>?xYjo+dhbMp~hsB?Pp#welkWq*Zcv{#AFUXxa+ z=_J8kSJ{Bd8+?ONLe*{tF8iA_qrF9R_O`S_y+aaghRPBYt9qAFLe*DQvG-_3d!OiR zrnExMA_+EIWgS(!MeE;js`de+gu3SQ(c?cKJ*{FN(v0>I(b>n+3iSy|uuoN%pjg#s zj1ua|jr#Y#t?l$C`N>vINztYTqzQs3#6xai=@{8x{MOX0-2! z&gMuf)b}L8eo)zf%J0D+86_0ExkP6_ODohol3>56EJ3lsUl}D-ZG%&q zY(C9s3y98slUAtTNrL^MvINBj|74UcK(ZIv}XJsue0XF0@Z>fSWA@+ zsC-|wVw6x9uhHwGHF{02SZkWm+7O+!l~${2M-pryl_e-P*q%{BP3Uv_oRb%RuVM?+ zjJ62T*`m@4wHQgT#Z{J|*x(Y35^75E3O-%2C22& zuq{+Jpz4q1s%^AXG7Xgd&{?I^8K zJCOw2S!D@|4fbS|P#u=k&pj2}g=Vx}iOzaSD^zciV0~29QU5sTcVm=L`Iy+9X0*OU zXM0F1)Se{4_EK4bVpWb&LKTPeJ{%K!(~Q=S=xiTph1!=S*nTP-Q2EKAKP^hAd^+z> zGui<}XaAB`r~^rY9i*}Z#i|Zwlu-FdJ%nboLy67?NGsGpl3;^W)={ng_H_x=U`7eW z1CJjL`XMx<9Y%C^xU@nIB?&f6C8(D9$&jDzM=(mLeA$nr8SN;dv!kUIYB))-5h_bi ztm+s>36(Fqry0$N&PGZrRFx#yD3uMUeA$m>lu-F}K8|Lz>|&J-sC+YA!YHBg>3k{8XqOS4T`sLqSC9m|Qe_E>4PM14q52g&ys=_e(~LHb z={^v2C^mQ@k%kC|312ql9|D?;iKe`i;+5O*ErDL3H+{v_efG3HFrA22|~`v#O^V zB~;7edBwGQhGw*9iO!yrR;cGmf=yLff?`!KFiNP|wN+K@MVirGB076nTA`+q1bam# zs8;#Wmsj;Fql79JrKjj?G^0%?I(uDOq23?~_NK}@>hFga8+?mVLe+l%>g|fXO*7g% zL}xRk73y7*VDG6cL9xO2870(=blGRpj5dqtY__yQeLxcILzNAvyuptcB~-peKc*S& z6QZ+Er4{Nkl3<^!EJ3lVFBm0M{bs1xmo%e&MRfMHv_gGD66{-*byV#Zt$h~dT7AbT zp~fuq?V9&r&hKB%p&9LaqO%{Q73xQlU_YrWL9wd2j1uaQ8%{mJzslF8pJ_&$M|Ad! zv_kz#5^TQ85>)G|Ent*T8^3wZPS3r>*QMWRM*E%U>e2A5-$P^~76Y&UM!jEXH!GujG7XDdo8R40;PE2%6& zvB8xYB~*{%73AyEDm0^YCOTVHTA@}W3AVaQP%VpRgI>M-jeZSA36;M)uSql7T101E zq!p?wNw97zOHkbHYcooyd^)c~GupaDXX{BTRCkhK>#Hn5ajiCBlu$D#9sM5vRBS_< z(KaGF+gMtmdXNO$L}ddiAM~3tN~rwZc{7^PHYYmULRz7=Bnh^a$`TYC+?r8B<epc!pXqO-lE70QtW+goJ`iVgN- zlu(n44Nj`qJ~X54OLVrMv_kbK3AVq=I_e(>{Q-;;Dxc2(q8aT#qO*gf73yG;V27wI zL9wbs86{LcCI--qHjwCSkhDS#CJ8n~WdkZd84hEVQ2BH|oMyD4L}$aK73v6*U`MJf zL9wc%7$sDz;&|qh%+WNX4JSGqA+1oykOcE8>!{X$`?>_m7$sC~gBcf?`!CF-oX>+5b&5+Q~#` zr${T*sU*QpQ`vyZm;H1`36)RhGiXMOL}z2A73xfqU}vc;L9wc{86{M{?Ej$|?Hr=B zbEOsPJd$ANtE{8`anN7DD53J{d?C$f7ZIIZEUi$NkOaF_WeJK^UB)P(@^|OUX-2z( z=3bas!lLfuOe>^_wZsC-}D&nTht>HGlAXp@M}9+Xz7 zhe(1wtg-~f1|MORQ0JVfAHD)SN;BGIqO-@O73y)4U`;AZP;BrCMhP{2i?xm(b}An~ zo}?LV3enk9(hBu7Nw8;B)={-5L+ur0RnIa?sH=<1&c~1EXhwUU=xnOALcKr|>_wF& zC|30nqlD@+^u2G!9Kh$RmuW_uMs)Uyv_ic~66`gV4XD~>XI0Y~B~*`%_49GXUZ)xD z4WhF*r4{Nel3;JEEJ3lVcNitq?zL4_YzEC}?-HH8C#_KLlLVWo5>)H_=*z2`#VDcj z*QMDsqkTYh_Mx;weMA!MW0iH(-w!V~_z9zgx+-1vPiaQ`jOgrhX@&ZNB-ocKOHgd^ zD@F;GFZ-m9-r)C)5^9rl*?*uJ?MI@spQIIPE=jPTRhFPw z)jUQCl`s1*G^71WbT(gFp%#z?`%PsXRl7xNpGCP=zcWgx>T5eMdr22oqdA-S6(!_s7#D68NP;E%?U#TovRhz1{W0X)^ zj$OZH-vxYKT8Jk8D1F<{8uVV zP;785MhW%Tbwl^~X=$EUU1;LJBI3W2R;aZ}@L#DcL9xMg870(K#Vg3yrS)jyzarwl zl2)h1mXmY}%XH)oVk`E=fb zCjKiT{wryP+L{FamC6zn*J@iv3AO(9kLi1A#VR!MUlH+NNh{P2B>1mXHlXrBzZ0W` z%HN%LriuTGi2q7jp>`#~f2Fbn#RhvbN~nA~_o0dZiirP8TA})q;J;E?f?|VvGD@g% zb9-Gow-;}my=daUBI3W2R;Yd?_^(tpplY`$-*NV3lu-G*^L{k(UlH+NNh{O=B>1mX zmY~?+fs7I=pUwx-#D7J^eIS_6aN(v|CO{t z9Y%uxN@X4OkAr?FqlC)G#4wupuZZ}sq!sEY68u*xOHiz8IHQCrP6?gPBWU8kBI3W2 zRwyIEf2Fbkm7fe%MhTUViBUB1UlH+NNh{RxB>1mXmY`VG35*h|_5qS7!-+KUUlH+N zNh{PzB>1mXf@+=L&iR3SGNXjb*Xk6S_^*igucQ^~bQ1hmDoaqTDl$r_{7GgkP5f6x z{8!Qnbv6n9E0qnXe67x5lu-Fvol6t{6%qfHv_f4#g8xco35r!+#3-Tir@o76;=dx| zzmis{%SiBFsjQ>^aXMeYD53I^dL>Q#S48|*(h4<>1pk%F5)`W%&nTht9|>Jc6aN(v z|CO{tO(4O4rLqLYv;77}2{mr=^g+vQ$akC@Y2v>k;=htss9Q+zU#VTNXfUlH+NNh{QyB>1mXmY~?+-HZ}yvEmPlFRa)-G^5>1#D68NQ1_GIzfxI(VuOVw66)ARXZ4;wm`{CEY2v>k;=hts zsFz6aU#SGuCcm9)8)Q|}7$sEg@0j4NeT63eDmLMEqCM3N?!a|CP!H zRNmkRj1sD5O2@ePAx->OMEqCM3iSyI{wtLwC|30uqlDU`wyKJKP80tX5&xC6LVZPo z|4JpO*2RVY-w#EcgWoVpsQd}}TblTAvINBje`S)Ow6Z=|B8tJN?M_QC&7QEvH_K!41Y38sQd}} zFPiwTS_{>jSfE;v;J;E?f?|WM7$sDGJGZ8Z|B8tJN?NU|9SQy`l_e-P*q%{Bow)e6 zzkYBUABq;HiT{d-|4Le+79+ubrLqB)pA1VdN~rt^c}bf1uZZ}sq!p?o3H~dUB`7wy z45Ng~Z|7xc;=dx|zmis{6-e-3sVqUU!A^`4ssmpJcV70fimgNw{}mDcm9#>2Cc%HD z5>%Ufi`HI2R<#5|B8tJN?M_IBEf&95|low)lm5r?8zvh@&47+{8uVVP^_vy zqlDTrJ&^aOiT{d-|4Le+4kW>Ur4p3R!P>Q|eM{xCAIvDBrgdI)I{#Gc5SsX}i1@Fh z6>1;}{wtLwC{{I?Q9}LpY{&KPy_Vld8A22P6%qfHv_cIf!GEQ)1l6``M=(mLdADzL z$DSRrCml%>{}mDcm9#<)C&7QEvH_KEhGQ5d)Zz0l8!`Vf9(|rB{wpH>D`|zQlHk8m zS%PAN$1+N&`NanBYqH~L;=dx|zmis{6G-r1sVqUU!7+>y%C~%}<>>Jo^noV+DsbZ(m#D7J^ewEj!{D8_u%z3@m~?~Ur8&}4J7!lRF|BI3W2R;Y<2 z_^(ug(v|wh+xd1z36r-n ze?`Q9C9P0XNbp~&EJ3lsrx_*GGI79ialo1QuV~`GBI3W2R;cGm@L#DcL9xLX7$wx< z<7~nKH{6H+iYERmBK|9Bg_=f!|4L;8D&JSHGD@iY9(;`^{wpH>D`|y#g9QJT$`TYC ze2Y;+RbKgN>sNl_o7dYk@m~?~Ur8&}yCnFpRFd?i%?6`Mxj)|Ex@m~?~ zUr8&}2PF8fRD#lHhQGg^S=C335^8obYAMuPuJWeJK^eZeT9_J3x@ z?h9_?Df%T%{8vQ$SJDdg4GI1$l?|xcWoK32F-oY*iofi@lVJ``{8vQ$SJDdgBMJU1 zl_e-vHJ4FB9Z*|U#eSxV|B8tJN?M_QCBc8C5>%VwYUxhTr}F|v36;;m-)Q2$BI3W2 zR;WKo@L#DcL9xMRZT=UP--FF*;=dx|zakc>RwVeZRFinhRWr`I<{kd2SX-L- zuZZ}sq}8h0liO;j!{AlUU1_f3vT2y!}2uoUlH+NNh?$*68u*x8&LV}yfUMN z%IDxJH1S^%@n1kS6{sBK|9Bh3Y|q z|4L;EidAjOD4{M*@4?Mz;=dx|zmis{ElKcSscb;yxAWGF5-Pt3x1ovuiirP8TA{Wh z!GEQ)1jVX$V3bh#d9@=={8vQ$SJDdAlLY^jN>J_cwW_TO)UJ#Y3a2DLSZ;dJ#D7J^ Qe?=TMfUh0=dtmqf0Z}4br2qf` literal 0 HcmV?d00001 diff --git a/examples/Su2ExactDiagonalization/Su2Engine.py b/examples/Su2ExactDiagonalization/Su2Engine.py new file mode 100644 index 000000000..021481ec3 --- /dev/null +++ b/examples/Su2ExactDiagonalization/Su2Engine.py @@ -0,0 +1,263 @@ +# Copyright 2020 The TensorNetwork Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Basic SU(2) representation theory data""" + +from math import floor +from math import factorial +from math import sqrt +from collections import defaultdict +import numpy as np +from tensornetwork.ncon_interface import ncon +from Su2Tensor import Su2Tensor +from Su2Tensor import FusionTree + +def myisinteger( + num: int) -> bool: + """ + Checks if num is an integer + """ + + val = 1 if num == floor(num) else 0 + return val + +####################################################################################### +def compatible( + j1: float, j2: float, j12: float) -> bool: + """ + Checks if total spins j1,j2,j12 are compatible with SU(2) fusion rules + """ + + val = True if j12 in np.arange(abs(j1-j2),j1+j2+1,1).tolist() else False + return val + +####################################################################################### +def wigner3j( + j1: float, m1: float, j2: float, m2: float, j12: float, m12: float) -> float: + """ + Computes the Wigner 3j symbol. j's at the total spin values and m's are the + corresponding spin projection values + """ + # compute Wigner 3j symbol + if not myisinteger(2*j1) or not myisinteger(2*j2) or not myisinteger(2*j12) \ + or not myisinteger(2*m1) or not myisinteger(2*m2) or not myisinteger(2*m12): + raise ValueError('Input values must be (semi-)integers!') + if not myisinteger(2*(j1-m1)): + raise ValueError('j1-m1 must be an integer!') + if not myisinteger(2*(j2-m2)): + raise ValueError('j2-m2 must be an integer!') + if not myisinteger(2*(j12-m12)): + raise ValueError('j12-m12 must be an integer!') + if not compatible(j1,j2,j12): + raise ValueError('The spins j1,j2,j12 are not compatible with the fusion rules!') + if abs(m1) > j1: + raise ValueError('m1 is out of bounds.') + if abs(m2) > j2: + raise ValueError('m2 is out of bounds.') + if abs(m12) > j12: + raise ValueError('m12 is out of bounds.') + if m1 + m2 + m12 != 0: + return 0 # m's are not compatible. So coefficient must be 0. + t1 = j2 - m1 - j12 + t2 = j1 + m2 - j12 + t3 = j1 + j2 - j12 + t4 = j1 - m1 + t5 = j2 + m2 + + tmin = max([0, t1, t2]) + tmax = min([t3, t4, t5]) + + wigner = 0 + for t in np.arange(tmin, tmax+1, 1).tolist(): + wigner = wigner + (-1)**t / (factorial(t) * factorial(t - t1) * factorial(t - t2)\ + * factorial(t3 - t) * factorial(t4 - t) * factorial(t5 - t)) + + wigner = wigner * (-1)**(j1-j2-m12) * sqrt(factorial(j1+j2-j12) * factorial(j1-j2+j12)\ + * factorial(-j1+j2+j12) / factorial(j1+j2+j12+1) * factorial(j1+m1)\ + * factorial(j1-m1) * factorial(j2+m2) * factorial(j2-m2) * factorial(j12+m12)\ + * factorial(j12-m12)) + + return wigner + +####################################################################################### +def clebschgordan( + j1: float, m1: float, j2: float, m2: float, j12: float, m12: float) -> float: + """ + Computes a Clebsch-Gordan coefficient. j's at the total spin values and m's are the corresponding + spin projection values + """ + return (-1)**(j1-j2+m12) * sqrt(2*j12+1) * wigner3j(j1,m1,j2,m2,j12,-m12) # Clebsch-Gordan coefficient from Wigner 3j + +####################################################################################### +def createCGtensor( + j1: float, j2: float, j12: float) -> np.array: + """ + Computes a Clebsch-Gordan tensor for the fusion j1xj2 -> j12 + """ + + C = np.zeros((floor(2*j1+1), floor(2*j2+1), floor(2*j12+1))) + if not compatible(j1,j2,j12): + return C + for m1 in np.arange(-j1, j1+1, 1).tolist(): + for m2 in np.arange(-j2, j2+1, 1).tolist(): + for m12 in np.arange(-j12, j12+1, 1).tolist(): + C[floor(j1+m1)][floor(j2+m2)][floor(j12+m12)] = clebschgordan(j1,m1,j2,m2,j12,m12) if m1+m2 == m12 else 0 + return C + +####################################################################################### +def rsymbol( + j1: float, j2: float, j12: float) -> float: + """ + Returns a R-symbol of SU(2) + """ + + C1 = createCGtensor(j1,j2,j12) + C2 = createCGtensor(j2,j1,j12) + C = ncon([C1, C2],[[1, 2, -1],[2, 1, -2]]) + return C[0][0] + +####################################################################################### +def racahDenominator( + t: float, j1: float, j2: float, j3: float, J1: float, J2: float, J3: float) -> float: + """ + Helper function used in wigner6j(). + """ + + return factorial(t-j1-j2-j3) * factorial(t-j1-J2-J3) * factorial(t-J1-j2-J3) * factorial(t-J1-J2-j3) \ + * factorial(j1+j2+J1+J2-t) * factorial(j2+j3+J2+J3-t) * factorial(j3+j1+J3+J1-t) + +####################################################################################### +def triangleFactor( + a: int, b: int, c: int) -> float: + """ + Helper function used in wigner6j(). + """ + + return factorial(a+b-c) * factorial(a-b+c) * factorial(-a+b+c)/factorial(a+b+c+1) + +####################################################################################### +def wigner6j( + j1: float, j2: float, j3: float, J1: float, J2: float, J3: float) -> float: + """ + Returns a Wigner 6j-symbol of SU(2). + """ + + tri1 = triangleFactor(j1, j2, j3) + tri2 = triangleFactor(j1, J2, J3) + tri3 = triangleFactor(J1, j2, J3) + tri4 = triangleFactor(J1, J2, j3) + + # Calculate the summation range in Racah formula + tmin = max([j1+j2+j3, j1+J2+J3, J1+j2+J3, J1+J2+j3]) + tmax = min([j1+j2+J1+J2, j2+j3+J2+J3, j3+j1+J3+J1]) + + Wigner = 0 + + for t in np.arange(tmin, tmax+1, 1).tolist(): + Wigner = Wigner + sqrt(tri1*tri2*tri3*tri4) * ((-1)**t) * factorial(t+1)/racahDenominator(t,j1,j2,j3,J1,J2,J3) + + return Wigner + +####################################################################################### +def fsymbol( + j1: float, j2: float, j3: float, j12: float, j23: float, j: float) -> float: + """ + Returns he recouping coefficient of SU(2), related to the Wigner 6j symbol by an overall factor. + """ + + if not compatible(j1,j2,j12) or not compatible(j2,j3,j23) or not compatible(j1,j23,j) or not compatible(j12,j3,j): + return 0 + + return wigner6j(j1,j2,j12,j3,j,j23) * sqrt(2*j12+1) * sqrt(2*j23+1) * (-1)**(j1+j2+j3+j) + +####################################################################################### +def isvalidindex( + index: dict) -> None: + """ + Checks for a valid SU(2) index, and raises an exception if its not. + index: a Dictionary of the type {total spin value: degeneracy dimension}. + """ + + if not isinstance(index, dict): + raise ValueError('SU(2) index should be a dictionary of the type: {total spin, degeneracy dimension}') + + for key in index: + if not myisinteger(2*key): + raise ValueError('Spin values in an index should be (semi-)integers.') + if not myisinteger(index[key]): + raise ValueError('Degeneracy values in an index should be integers.') + +######################################################################################### +def fuseindices( + ia: dict, ib: dict) -> dict: + """ + Fuses two SU(2) indices ia, ib according to the fusion rules. + ia, ib: Dictionaries of the type {total spin value: degeneracy dimension}. + """ + + if len(ia) == 0 or len(ib) == 0: + raise ValueError('Input indices cannot be empty.') + isvalidindex(ia) + isvalidindex(ib) + + iab = defaultdict(int) # initialize fused index to empty default dictionary + for ja in ia: + da = ia[ja] + for jb in ib: + db = ib[jb] + for jjab in np.arange(abs(ja-jb), ja+jb+1, 1).tolist(): + if len(iab) == 0: + iab[jjab] = da*db + else: + iab[jjab] += da*db + + return iab + +######################################################################################## +def fuse( + ia: dict, ib: dict, totalSectors: list=[]) -> Su2Tensor: + """ + Creates a 3-index tensor X that two SU(2) indices according to the fusion rules. + ia, ib: Dictionaries of the type {total spin value: degeneracy dimension}. + totalSectors (optional): list of total sectors to keep. [total spin values] + """ + iab = fuseindices(ia, ib) + + # remove some total sectors if targetSectors are specified + if totalSectors: + iabCopy = dict(iab) # shallow copy should suffice + for jab in iab: + if not jab in totalSectors: + iabCopy.pop(jab) # remove unwanted total sectors + iab = iabCopy # update iab + + dabGone = defaultdict(int) # dictionary index by jab to keep track to degeneracies as they are filled in X + + X = Su2Tensor([iab, ia, ib]) + for ja in ia: + da = ia[ja] + for jb in ib: + db = ib[jb] + for jab in iab: + if not compatible(ja,jb,jab): + continue + dab = iab[jab] + temp = np.zeros((da*db,dab)) + temp[0:db*da, dabGone[jab]:db*da+dabGone[jab]] = np.identity(db*da) + temp = temp.reshape((db, da, dab), order='F') # Specifying order=F means a column-major reshape, similar to MATLAB. Default in python is row-major + temp = temp.transpose(2,1,0) + X.addblock(FusionTree((jab, ja, jb)),temp) + dabGone[jab] += da*db + + return X \ No newline at end of file diff --git a/examples/Su2ExactDiagonalization/Su2ExactDiagonalization.py b/examples/Su2ExactDiagonalization/Su2ExactDiagonalization.py new file mode 100644 index 000000000..fbcf9f9e6 --- /dev/null +++ b/examples/Su2ExactDiagonalization/Su2ExactDiagonalization.py @@ -0,0 +1,277 @@ +# Copyright 2020 The TensorNetwork Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" Performs exact diagonalization of SU(2)-symmetric Hamiltonians by exploiting the symmetry. + - Diagonalization can be restricted to specified total spin sectors. + - Uses pre-computation of contraction factors resulting from contracting Clebsch-Gordan tensors; + factors are computed on first run of the code and stored in a file. + - Provides a preset of popular SU(2)-symmetric Hamiltonian densities + - Assumptions: 1) Hamiltonian density must be a 2-site operator. 2) Presently only implements open + boundary condition. 3) Only uses ncon() facility of the TensorNetwork library.""" + +import numpy as np +from numpy import linalg as LA +from math import floor +from scipy.sparse.linalg import eigs + +from tensornetwork.ncon_interface import ncon + +from Su2Tensor import Su2Tensor +from Su2Tensor import FusionTree +import Su2Engine + +# globals +total, left, right = (0, 1, 2) # some constants (index values) + +######################################################################################################################## +def Su2ExactDiagonalization( + N: 'int', d: 'dict', h: Su2Tensor, targetSectors: list=[], + tableFileName: 'string' = 'FT_Su2CoarseGrain2To1Site.npy', + doCreateTable: 'bool' =False) -> tuple[dict, dict, np.array]: + """ + Performs exact diagonalization in specified targetSectors. Returns the low-energy states and their energies in + these sectors. + N: number of lattice sites. + d: SU(2) index describing each lattice site. Must be a dictionary of the type {total spin: degeneracy dimension} + h: 2-site Hamiltonian density. Must be an object of SU2Tensor class + targetSectors: Specifies which total spin sectors (of the entire lattice) to diagonalize and how many states + required in each sector. Must be dictionary of the type: {total spin, number of states} + If None diagonalizes all the sectors + doCreateTable: set to True if the factor tables required by this function have not been created yet. + tableFileName: file to store the table when table creation is demanded. + """ + + # check validity of inputs + if not Su2Engine.myisinteger(N): + raise ValueError('Lattice size must be a natural number!') + if N == 0 or N == 1: + raise ValueError('Lattice size must be greater than equal to 2!') + Su2Engine.isvalidindex(d) + if not isinstance(h,Su2Tensor): + raise ValueError('Hamiltonian density must be a SU(2) tensor.') + + # create blocking isometries from bottom (physical index) to top(total index) + w = [] # list of blocking isometries + bond = d + for k in range(N-2): + X = Su2Engine.fuse(bond, d) + w.append(X) + bond = X.getIndices()[0] # set bond to new total index + + # Create the final blocking isometry specially + if targetSectors: # the user has asked to truncate final index + X = Su2Engine.fuse(bond, d, list(targetSectors.keys())) # truncate total spin sectors + else: + X = Su2Engine.fuse(bond, d) # otherwise the full blocking tensor + w.append(X) + bond = X.getIndices()[0] # bond is now the total spin basis on the entire lattice + + if doCreateTable: + factorTable = createFactorTable_Su2CoarseGrain_2To1Site([[w[N-4],w[N-3]], [w[N-3],w[N-2]]],w[0], tableFileName) + else: + factorTable = np.load(tableFileName, allow_pickle='TRUE').item() + + su2H = Su2Tensor([bond, bond]) # Total Hamiltonian. + wd = w[0] + for term in range(N-1): + # first coarse-grain local term to a 1-site operator. Fist term has to + # be treated differently from the others. + + if term == 0: + hh = h # already blocked to a one-site operator + else: + hh = Su2CoarseGrain_2To1site(w[term-1], w[term], wd, h, factorTable) + + # now coarse-grain this one-site operator through the remaining linear fusion tree + for m in range(term+1,N-1): + hh = Su2CoarseGrain_1To1site(w[m],hh) + + # hh is now the one-site version of the local hamiltonian term. Add to current sum + hhBlocks = hh.getAllBlocks() + for fusionTree in hhBlocks: + su2H.addblock(fusionTree,hhBlocks[fusionTree]) + + # Now let's diagonalize + spectrum = dict() # dictionary of the type: {total spin: spectra} + eigenvectors = dict() # dictionary of the type: {total spin: eigenvectors} + fullspectrum = None # full spectrum (including degneracies) as a NumPy array + su2HBlocks = su2H.getAllBlocks() + for fusionTree in su2HBlocks: + degdim = 2*fusionTree[0]+1 + if targetSectors: # diagonalize only the sectors that the user wants + temp, ev = eigs(su2HBlocks[fusionTree], k=targetSectors[fusionTree[0]], which='SR') + else: # do full digonalization of all the blocks + temp = LA.eigvals(su2HBlocks[fusionTree]) # full diagonalization + ev = None + spectrum[(fusionTree[0])] = temp + eigenvectors[(fusionTree[0])] = ev + temp = np.kron(temp,np.ones((1,floor(degdim)))) + if fullspectrum is None: + fullspectrum = temp + else: + fullspectrum = np.append(fullspectrum, temp) + fullspectrum.sort() + + return spectrum, eigenvectors, fullspectrum + +######################################################################################################################## +def Su2CoarseGrain_2To1site( + w1: Su2Tensor, w2: Su2Tensor, wd: Su2Tensor, h: Su2Tensor, factorTable: dict) -> Su2Tensor: + """ + Coarse-grains 2-site Ham h to a 1-site term through isometries w1,w2 inside a linear fusion tree + w2 is the blocking isometry from the higher level. h,hL,w1,w2 are SU(2)-symmetric tensors. + """ + + hh = Su2Tensor(h.getIndices()) # the coarse-grained hamiltonian. Even though h is a 4-index tensor it is + # stored as a 2-index matrix. So indices of h are the blocked 2-site indices + w2Blocks = w2.getAllBlocks() + w1Blocks = w1.getAllBlocks() + wdBlocks = wd.getAllBlocks() + hBlocks = h.getAllBlocks() + for ft in factorTable: # iterate through all fusion trees in table + c1 = FusionTree((ft[0], ft[1], ft[4])) + c2 = FusionTree((ft[1], ft[2], ft[3])) + c3 = FusionTree((ft[5], ft[3], ft[4])) + c4 = FusionTree((ft[5], ft[6], ft[7])) + c5 = FusionTree((ft[8], ft[2], ft[6])) + c6 = FusionTree((ft[0], ft[8], ft[7])) + ch = FusionTree((ft[5], ft[5])) + + if not c1 in w2Blocks or not c2 in w1Blocks or not c3 in wdBlocks or not c4 in wdBlocks or not c5 in w1Blocks \ + or not c6 in w2Blocks or not ch in hBlocks: + continue + temp = ncon([w2Blocks[c1],w1Blocks[c2],wdBlocks[c3],hBlocks[ch],wdBlocks[c4],w1Blocks[c5],w2Blocks[c6]],\ + [[-1,2,4],[2,8,3],[1,3,4],[1,9],[9,6,7],[5,8,6],[-2,5,7]]) + hh.addblock(FusionTree((ft[0],ft[0])),factorTable[ft]*temp) + + return hh + +######################################################################################################################## +def Su2CoarseGrain_1To1site( + w: Su2Tensor, h: Su2Tensor) -> Su2Tensor: + """ + Coarse-grains (blocks) 1-site operator h through isometry w. + w: the blocking isometry + h: the 1-site SU(2)-symmetric operator + """ + + hh = Su2Tensor(h.getIndices()) # the coarse-grained hamiltonian. Even though h is a 4-index tensor it is + # stored as a 2-index matrix. So indices of h are the blocked 2-site indices + wBlocks = w.getAllBlocks() # a dictionary {FusionTree:block} + hBlocks = h.getAllBlocks() # a dictionary {FusionTree:block} + for ft in wBlocks: + ch = FusionTree((ft[1],ft[1])) + if not ch in hBlocks: + continue + temp = ncon([wBlocks[ft],hBlocks[ch],wBlocks[ft]],[[-1, 1, 3],[1, 2],[-2, 2, 3]]) + hh.addblock(FusionTree((ft[0],ft[0])),temp) + + return hh + +######################################################################################################################## +def createFactorTable_Su2CoarseGrain_2To1Site( + W: 'list of pairs of Su2Tensors', wd: Su2Tensor, tableFileName: str) -> dict: + """ + (Precomputation.) Creates the factors and compatible chargesectors for the function Su2CoarseGrain_2To1Site() + W: a list of pairs of SU(2) tensors + wd: a SU(2) tensor (the isometry that blocks 2 lattice sites) + tableFileName: file name to store the created factor table + """ + + FactorTable = dict() # Dictionary from fusion trees to factors + for k in range(len(W)): + w1 = W[k][0] + w2 = W[k][1] + for c1 in w2.getAllBlocks(): + for c2 in w1.getAllBlocks(): + if c2[total] != c1[left]: + continue + for c3 in wd.getAllBlocks(): + if c2[right] != c3[left] or c1[right] != c3[right] or not Su2Engine.compatible(c2[left], c3[total],\ + c1[total]): + continue + for c4 in wd.getAllBlocks(): + if c4[total] != c3[total]: + continue + for c5 in w1.getAllBlocks(): + if c5[right] != c4[left] or c5[left] != c2[left] or not Su2Engine.compatible(c5[total],\ + c4[right],c1[total]): + continue + for c6 in w2.getAllBlocks(): + if c6[left] != c5[total] or c6[right] != c4[right] or c6[total] != c1[total]: + continue + + C1 = Su2Engine.createCGtensor(c1[left], c1[right], c1[total]) + C2 = Su2Engine.createCGtensor(c2[left], c2[right], c2[total]) + C3 = Su2Engine.createCGtensor(c3[left], c3[right], c3[total]) + C4 = Su2Engine.createCGtensor(c4[left], c4[right], c4[total]) + C5 = Su2Engine.createCGtensor(c5[left], c5[right], c5[total]) + C6 = Su2Engine.createCGtensor(c6[left], c6[right], c6[total]) + Cmat = ncon([C1,C2,C3,C4,C5,C6],[[1,3,-1],[7,2,1],[2,3,8],[5,6,8],[7,5,4],[4,6,-2]]) + FactorTable[FusionTree((c1[total],c1[left],c2[left],c2[right],c1[right],c3[total],\ + c4[left],c4[right],c5[total]))] = Cmat[0][0] + + np.save(tableFileName, FactorTable) + return FactorTable + +######################################################################################################################## +def Su2Create2SiteHamiltonianDensity( + whichHam: str) -> Su2Tensor: + """ + Creates 2-site Hamiltonian density as a SU(2) tensor. d is the SU(2) index for each lattice site. + whichHam is a string to select from a preset of SU(2)-symmetry Hamiltonian densities: + 'afmHeisenberg', 'spin1Heisenberg', 'aklt'""" + if whichHam == 'afmHeisenberg': + d2 = {0:1,1:1} + h = Su2Tensor([d2, d2]) + block0 = np.ones((1,1)) + block0[0][0] = -3 + block1 = np.ones((1,1)) + h.addblock(FusionTree((0, 0)), block0) + h.addblock(FusionTree((1, 1)), block1) + elif whichHam == 'identity': # for testing purposes + d2 = {0: 1, 1: 1} + h = Su2Tensor([d2, d2]) + block0 = np.ones((1, 1)) + block1 = np.ones((1, 1)) + h.addblock(FusionTree((0, 0)), block0) + h.addblock(FusionTree((1, 1)), block1) + elif whichHam == 'spin1Heisenberg': + d2 = {0: 1, 1: 1, 2: 1} + h = Su2Tensor([d2, d2]) + block0 = np.ones((1, 1)) + block0[0][0] = -2 + block1 = np.ones((1, 1)) + block1[0][0] = -1 + block2 = np.ones((1, 1)) + block2[0][0] = 1 + h.addblock(FusionTree((0, 0)), block0) + h.addblock(FusionTree((1, 1)), block1) + h.addblock(FusionTree((2, 2)), block2) + elif whichHam == 'aklt': + d2 = {0: 1, 1: 1, 2: 1} + h = Su2Tensor([d2, d2]) + block0 = np.ones((1, 1)) + block0[0][0] = -2/3 + block1 = np.ones((1, 1)) + block1[0][0] = -2/3 + block2 = np.ones((1, 1)) + block2[0][0] = 4/3 + h.addblock(FusionTree((0, 0)), block0) + h.addblock(FusionTree((1, 1)), block1) + h.addblock(FusionTree((2, 2)), block2) + elif whichHam == 'blbqHeisenberg': + pass + + return h diff --git a/examples/Su2ExactDiagonalization/Su2Tensor.py b/examples/Su2ExactDiagonalization/Su2Tensor.py new file mode 100644 index 000000000..f38ad62bb --- /dev/null +++ b/examples/Su2ExactDiagonalization/Su2Tensor.py @@ -0,0 +1,90 @@ +# Copyright 2020 The TensorNetwork Authors +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Definition of classes FusionTree and Su2Tensor. (Very basic at present, only handles 2- and 3- index tensors)""" +import numpy as np + +class FusionTree: + """ + Class to encapsulate a FusionTree. Currently used only for 2-index and 3-index Su2Tensors. In these cases FusionTree + is essentially a wrapper on a list of charges (2 or 3 charges) + """ + def __init__(self, charges): + self.charges = charges # charges is simply a tuple for now. (Should be replaced by a decorated tree in future) + def __hash__(self): + """ + To use objects of FusionTree as a key in a dictionary + """ + return hash(self.charges) + + def __eq__(self, other): + """ + To use objects of FusionTree as a key in a dictionary. + """ + return self.charges == other.charges + + def __ne__(self, other): + """ + To use objects of FusionTree as a key in a dictionary. Not strictly necessary, but to avoid having + both x==y and x!=y True at the same time. + """ + return not(self == other) + + def __getitem__(self, key): + """ + Override the built-in [] operator, so we can access elements of charges directly using the object + Currently, the class FusionTree is a wrapper to a list charges. + """ + return self.charges[key] + +class Su2Tensor: + """ + Only for 2- and 3-index SU(2) tensors for now! Currently essentially a wrapper for a dict: {FusionTree,block} + """ + def __init__( + self, indices: list): + self.blocks = dict() # blocks is a list indexed by keys which are FusionTrees + self.indices = indices + + def addblock( + self, fusiontree: FusionTree, block: np.array): + """ + Adds a block indexed by the given fusiontree. + """ + if fusiontree in self.blocks: + self.blocks[fusiontree] += block + else: + self.blocks[fusiontree] = block + + def getblock( + self, fusiontree: FusionTree): + """ + Retrieves a block indexed by the given fusiontree. + """ + if fusiontree in self.blocks: + return self.blocks[fusiontree] + else: + raise KeyError('Fusion tree not found in the tensor.') + + def getAllBlocks(self): + """ + Retrieves the whole block dictionary. + """ + return self.blocks + + def getIndices(self): + """ + Retrieves all the indices of the tensor. + """ + return self.indices diff --git a/examples/Su2ExactDiagonalization/main.py b/examples/Su2ExactDiagonalization/main.py new file mode 100644 index 000000000..b3bd7bfb9 --- /dev/null +++ b/examples/Su2ExactDiagonalization/main.py @@ -0,0 +1,113 @@ +""" +Test cases for SU(2) Exact Diagonalization +""" +import numpy as np +from tensornetwork.ncon_interface import ncon +import Su2Engine +from Su2Tensor import FusionTree +import Su2ExactDiagonalization as SD +from scipy.sparse.linalg import eigs + +######################################################################################################################## +if __name__ == '__main__': + + ###### Test1 : diagonalize spin-1 Heisenberg model on a lattice of N=4 sites ####################################### + print('Test1: spin-1 Heisenberg model on a lattice of N=4 sites') + N = 4 # number of lattice sites + d = {1:1} # SU(2) charge on each lattice site = spin 1 + h = SD.Su2Create2SiteHamiltonianDensity('spin1Heisenberg') + + # full diagonalization. When running Su2ExactDiagonalization() for the first time for given N,d set doCreateTable + # to True in order to generate the factor tables for use in subsequent calls (even for different ham, see Test2). + # The table for N also work for any small M < N. + spec = SD.Su2ExactDiagonalization(N,d,h,tableFileName = 'spinOneFactorTable.npy',doCreateTable = True) + SD.printSpectrum(spec[0]) # prints sector-wise spectrum, but degeneracies are suppressed + print('===========================================================') + ###### End of Test1 : diagonalize spin-1 Heisenberg model on a lattice of N=4 sites ################################ + + ###### Test2 : diagonalize spin-1 AKLT model on a lattice of N=4 sites ############################################# + print('Test2: spin-1 AKLT model on a lattice of N=4 sites') + N = 4 # number of lattice sites + d = {1:1} # SU(2) charge on each lattice site = spin 1 + h = SD.Su2Create2SiteHamiltonianDensity('aklt') + + # full diagonalization. Reuse factor table from Test1 + spec = SD.Su2ExactDiagonalization(N,d,h,tableFileName = 'spinOneFactorTable.npy',doCreateTable = False) + SD.printSpectrum(spec[0]) # prints sector-wise spectrum, but degeneracies are suppressed + + # diagonalize to find 2 lowest states in total spin = 0, 3 states in total spin = 1, 2 states in total spin 2 + spec = SD.Su2ExactDiagonalization(N, d, h, targetSectors = {0:2, 1:3, 2:2}, tableFileName='spinOneFactorTable.npy', + doCreateTable=False) # reuse factor tables from previous runs + SD.printSpectrum(spec[0]) # prints sector-wise spectrum, but degeneracies are suppressed + print('===========================================================') + ###### End of Test2 : diagonalize spin-1 AKLT model on a lattice of N=4 sites ###################################### + + ###### Test3 : diagonalize spin-1/2 AFM Heisenberg model on a lattice of N=8 sites ################################# + print('Test2: spin-1/2 AFM Heisenberg model on a lattice of N=8 sites') + N = 8 # number of lattice sites + d = {0.5: 1} # SU(2) charge on each lattice site = spin 1/2 + h = SD.Su2Create2SiteHamiltonianDensity('afmHeisenberg') + + # full diagonalization. Create factor tables for spin 1/2 chain with N=8 sites and store in specified file. + spec = SD.Su2ExactDiagonalization(N, d, h, tableFileName='spinHalfFactorTable.npy', doCreateTable=True) + SD.printSpectrum(spec[0]) # prints sector-wise spectrum, but degeneracies are suppressed + print('===========================================================') + ###### End of Test3 : diagonalize spin-1 AKLT model on a lattice of N=4 sites ###################################### + + ###### Test4 : Calculate Wigner and Clebsch-Gordan coefficients of SU(2) ########################################### + print('Demonstration of functions to calculate Wigner and Clebsch-Gordan coefficients of SU(2)') + w3j = Su2Engine.wigner3j(0.5, 0.5, 0.5, -0.5, 0, 0) + cg = Su2Engine.clebschgordan(0.5, 0.5, 0.5, -0.5, 0, 0) + cg1 = Su2Engine.clebschgordan(1, -1, 1, 0, 1, -1) + cg2 = Su2Engine.clebschgordan(1, -1, 1, 0, 2, -2) + C = Su2Engine.createCGtensor(0.5, 0.5, 1) + factor = Su2Engine.fsymbol(1, 1, 1, 1, 1, 1) + factor1 = Su2Engine.fsymbol(1, 2, 0.5, 3, 1.5, 2.5) + print('===========================================================') + ###### End of Test4 : Calculate Wigner and Clebsch-Gordan coefficients of SU(2) #################################### + + ###### Test5 : Check fusion rules of SU(2) ######################################################################### + print('Demonstration of functions to check fusion rules of SU(2)') + res = Su2Engine.compatible(j1=0.5, j2=0.5, j12=0.5) + res = Su2Engine.compatible(j1=1, j2=1, j12=1) + res = Su2Engine.compatible(j1=1, j2=0.5, j12=1) + print('===========================================================') + ###### End of Test5 : Check fusion rules of SU(2) ################################################################## + + ###### Test5 : Fusing two SU(2) indices ############################################################################ + print('Demonstration of fusing two SU(2) indices') + ia = {0.5:1} + ib = {0.5:1} + ic = {0.5:3, 2.5:1.3} + + try: + Su2Engine.isvalidindex(ia) + except: + print('Invalid index') + + try: + Su2Engine.isvalidindex(ic) + except: + print('Invalid index') + + iab = Su2Engine.fuseindices(ia,ib) # fuses indices ia and ib. index iab has the fused spin and degeneracy values + iiab = Su2Engine.fuseindices(iab,iab) # another example + + ia = {0:1, 1:1} + ib = {0:1, 1:1} + X = Su2Engine.fuse(ia,ib,[0,1]) # X is the fusion tensor. Its a Su2Tensor that can be contracted with another + # Su2Tensor inorder to fuse two indices and reshape that tensor + block000 = X.getblock(FusionTree((0, 0, 0))) # should be some assignments of 0's and 1's to give an injective map + # from input indices to fused index + block011 = X.getblock(FusionTree((0, 1, 1))) + block101 = X.getblock(FusionTree((1, 0, 1))) + block110 = X.getblock(FusionTree((1, 1, 0))) + block111 = X.getblock(FusionTree((1, 1, 1))) + + ia = {0.5:1, 1:2, 3:4} + ib = {0.5:1, 5:5, 2:1} + X = Su2Engine.fuse(ia,ib) + block0 = X.getblock(FusionTree((0, 0.5, 0.5))) + block1 = X.getblock(FusionTree((1, 1, 2))) + print('===========================================================') + ###### End of Test5 : Fusing two SU(2) indices ##################################################################### diff --git a/examples/Su2ExactDiagonalization/spinOneFactorTable.npy b/examples/Su2ExactDiagonalization/spinOneFactorTable.npy new file mode 100644 index 0000000000000000000000000000000000000000..cc16d159de51358c1b0ba4ed74bef32e3ca3008d GIT binary patch literal 19569 zcmbuGcl_2<|Hs=UJJ}(#GP1JAWky9fLflr1GRk>(h3L4=6K;Etgl;o?WyP1wjF1^B zdy|nU>-T(puIrqi%j5I=>-W9gZ|~vt?3~Z_xjyH*>ea=4_t|~V{hK#Cw%LhY1`ioI z$hvIVxyzV6x~$u|%V8s|Y6FHJGQtKAQM>1Wqec#)KXT}RV}?+7?Y3^`(E0!O)Bfa^ zgNBbf`k2w(2aT{H-H#r1RCU+@vjL;qA7VoWjTk<%YNG~K+xup<|Em4)!L_Mx-ly3a zzQvBs2l|#{eJkf%ckR=>Tm3&8TpfMPkoLY!{}#=fHQQrz-}bCC$NF~8FVwY9%l^&j zjq&aK_U+sIFMf*9r&9F`4?OE2zewN|`a_9~LYBH~63!{EYE#Ac zioUn=eG1`j!S7xZ)~@E??xwy;m|w>|qTkc`y$YcVe(y3NOYN70F1;N4M8B`|`xV0e z!S7!tWT^)v;p5*AT>RX#-)r*!ivB?74=RKQ2Y*PTFy8=&CSiUZ2Sh*6`9XznaPULQ zge>*2B>Y`(%D=ef4v&7Q^TP_^5y2l>CS<8cCEB;ifF_5b38N!$PzMt_m>7Z<`yg1@v;SpRja{<0+e zv-lN1ulUQOzry(|3*lA4UtK2TR~(mwx7L2eioYiM@y=gc2(JtNdI=jYjNZysKOqT! zOBd(A(cj?wjfL>0;BQXC`s2BH0wBC434cw(Tce-o{B4Es_TcYG!rJcg9`8)T+8%j1 z?u!0y=kFTe=!Mb zS2h2NFC}5VI4?&(&G}af;j6*FmV`A){S~eH=}CC`mHL@l@vle!hVyS0!ncBdyG+PZ z-$}xE=e;qZUwpO6&xrnA=ie)Y?*~7#OxUvOXC>izo$p>~$hogq{OssIaQ?$W_)+j5 zHwyFH>XRg#`jUR$ulP@+|IGQ%3*i^Re_1AEsb3}G&Q_Uw(o|05ucQCQ`ELv1cfrpo z6SCCrlkoKo&OC3~elwf=57GbV{7;2&Zty=h3iHb`FA0Ynrk}Se{+H-~b$)&!ToC+k zWkQzvdlELgdv2$WU$e(QqW{zRza+F~e2NOzoLDArS+y375}|f^pQ>0(n$cPjowb%$ zs5T_Q+Nvx*#?+UdKTKIWMhTTZkk~>rqqQeGTUc7579k0?sLB#l^QtYzD528PW{cB| zwgl1HlF|y*fh5>cDoao;s@9QFLZy!owlvLX%MhI{E3Ht=kpx>_Wdo{qJ2k7?3XBpe zeT=acX-4ZrbheVTLaj^^Y!#IyC^p!cQ9`8y(N?7yZ8f5^)uk0`4U%AMsw_dV!L=AA zR66Ob3(aUMhR7X2K3`9wgt^-TN0gZC9P0flLXsF zWgV4n(&9C&+O~`mDxH{Cp&4yEqO0~PU^}WTL9wcx7$sCXLv3f8(Rvb{?INvE zyOIR!rLqB)-uu>@Q9`9t)%wtkwj0se?$Qd?mn7I8DoaqTYEMQ9Mas{s72As}nj<>f zn^>Uwkp$aEWgXSBIB;vL0<|xrgsPoNb1JqU&1m{==8P|+HPitl!TzNZREvL}s8u_V zQ9`AY(GH>+?O>v_L!{LP4+Xr~gLohGeNr;`LbLuCUh9m*CNB~&_IZ7j`b zXA+&AC9P0rlLY&Z$`TYCJcm(2?Okl}fr_0=GunAXXXi^R)CDBLE>u~9VuKekN~q_H z4Nk4t#WbT`LUeYiv_f4*66|u74XAuyUBM`!(xGfu(u{T$(b?6~3N?--*flCkP;78K zqlCJn_!jd~#jd3p?K+~f>!lTH0!gs{sw_dV!5bJQRP7aHgE!KQb`#Or&C&{W3rVnB zRn}3pCqwNOWK|OxCDgv@zPgQOwA+c!?vPffJ4u4wrLqLYs_tf#P=_!0b7oj)R+HUB zGupjGXZJ}f)cqvE9#GkU%FnAwj1sDUv8qoi_8`q@4-uU`EUi$FkOX^FWeJK^O=gr( z!)mLl*kd%KJx+AiB&|?SkOX^DC8(D9(U(^>g;7G)o($ht>?xYjo+dhbMp~hsB?Pp#welkWq*Zcv{#AFUXxa+ z=_J8kSJ{Bd8+?ONLe*{tF8iA_qrF9R_O`S_y+aaghRPBYt9qAFLe*DQvG-_3d!OiR zrnExMA_+EIWgS(!MeE;js`de+gu3SQ(c?cKJ*{FN(v0>I(b>n+3iSy|uuoN%pjg#s zj1ua|jr#Y#t?l$C`N>vINztYTqzQs3#6xai=@{8x{MOX0-2! z&gMuf)b}L8eo)zf%J0D+86_0ExkP6_ODohol3>56EJ3lsUl}D-ZG%&q zY(C9s3y98slUAtTNrL^MvINBj|74UcK(ZIv}XJsue0XF0@Z>fSWA@+ zsC-|wVw6x9uhHwGHF{02SZkWm+7O+!l~${2M-pryl_e-P*q%{BP3Uv_oRb%RuVM?+ zjJ62T*`m@4wHQgT#Z{J|*x(Y35^75E3O-%2C22& zuq{+Jpz4q1s%^AXG7Xgd&{?I^8K zJCOw2S!D@|4fbS|P#u=k&pj2}g=Vx}iOzaSD^zciV0~29QU5sTcVm=L`Iy+9X0*OU zXM0F1)Se{4_EK4bVpWb&LKTPeJ{%K!(~Q=S=xiTph1!=S*nTP-Q2EKAKP^hAd^+z> zGui<}XaAB`r~^rY9i*}Z#i|Zwlu-FdJ%nboLy67?NGsGpl3;^W)={ng_H_x=U`7eW z1CJjL`XMx<9Y%C^xU@nIB?&f6C8(D9$&jDzM=(mLeA$nr8SN;dv!kUIYB))-5h_bi ztm+s>36(Fqry0$N&PGZrRFx#yD3uMUeA$m>lu-F}K8|Lz>|&J-sC+YA!YHBg>3k{8XqOS4T`sLqSC9m|Qe_E>4PM14q52g&ys=_e(~LHb z={^v2C^mQ@k%kC|312ql9|D?;iKe`i;+5O*ErDL3H+{v_efG3HFrA22|~`v#O^V zB~;7edBwGQhGw*9iO!yrR;cGmf=yLff?`!KFiNP|wN+K@MVirGB076nTA`+q1bam# zs8;#Wmsj;Fql79JrKjj?G^0%?I(uDOq23?~_NK}@>hFga8+?mVLe+l%>g|fXO*7g% zL}xRk73y7*VDG6cL9xO2870(=blGRpj5dqtY__yQeLxcILzNAvyuptcB~-peKc*S& z6QZ+Er4{Nkl3<^!EJ3lVFBm0M{bs1xmo%e&MRfMHv_gGD66{-*byV#Zt$h~dT7AbT zp~fuq?V9&r&hKB%p&9LaqO%{Q73xQlU_YrWL9wd2j1uaQ8%{mJzslF8pJ_&$M|Ad! zv_kz#5^TQ85>)G|Ent*T8^3wZPS3r>*QMWRM*E%U>e2A5-$P^~76Y&UM!jEXH!GujG7XDdo8R40;PE2%6& zvB8xYB~*{%73AyEDm0^YCOTVHTA@}W3AVaQP%VpRgI>M-jeZSA36;M)uSql7T101E zq!p?wNw97zOHkbHYcooyd^)c~GupaDXX{BTRCkhK>#Hn5ajiCBlu$D#9sM5vRBS_< z(KaGF+gMtmdXNO$L}ddiAM~3tN~rwZc{7^PHYYmULRz7=Bnh^a$`TYC+?r8B<epc!pXqO-lE70QtW+goJ`iVgN- zlu(n44Nj`qJ~X54OLVrMv_kbK3AVq=I_e(>{Q-;;Dxc2(q8aT#qO*gf73yG;V27wI zL9wbs86{LcCI--qHjwCSkhDS#CJ8n~WdkZd84hEVQ2BH|oMyD4L}$aK73v6*U`MJf zL9wc%7$sDz;&|qh%+WNX4JSGqA+1oykOcE8>!{X$`?>_m7$sC~gBcf?`!CF-oX>+5b&5+Q~#` zr${T*sU*QpQ`vyZm;H1`36)RhGiXMOL}z2A73xfqU}vc;L9wc{86{M{?Ej$|?Hr=B zbEOsPJd$ANtE{8`anN7DD53J{d?C$f7ZIIZEUi$NkOaF_WeJK^UB)P(@^|OUX-2z( z=3bas!lLfuOe>^_wZsC-}D&nTht>HGlAXp@M}9+Xz7 zhe(1wtg-~f1|MORQ0JVfAHD)SN;BGIqO-@O73y)4U`;AZP;BrCMhP{2i?xm(b}An~ zo}?LV3enk9(hBu7Nw8;B)={-5L+ur0RnIa?sH=<1&c~1EXhwUU=xnOALcKr|>_wF& zC|30nqlD@+^u2G!9Kh$RmuW_uMs)Uyv_ic~66`gV4XD~>XI0Y~B~*`%_49GXUZ)xD z4WhF*r4{Nel3;JEEJ3lVcNitq?zL4_YzEC}?-HH8C#_KLlLVWo5>)H_=*z2`#VDcj z*QMDsqkTYh_Mx;weMA!MW0iH(-w!V~_z9zgx+-1vPiaQ`jOgrhX@&ZNB-ocKOHgd^ zD@F;GFZ-m9-r)C)5^9rl*?*uJ?MI@spQIIPE=jPTRhFPw z)jUQCl`s1*G^71WbT(gFp%#z?`%PsXRl7xNpGCP=zcWgx>T5eMdr22oqdA-S6(!_s7#D68NP;E%?U#TovRhz1{W0X)^ zj$OZH-vxYKT8Jk8D1F<{8uVV zP;785MhW%Tbwl^~X=$EUU1;LJBI3W2R;aZ}@L#DcL9xMg870(K#Vg3yrS)jyzarwl zl2)h1mXmY}%XH)oVk`E=fb zCjKiT{wryP+L{FamC6zn*J@iv3AO(9kLi1A#VR!MUlH+NNh{P2B>1mXHlXrBzZ0W` z%HN%LriuTGi2q7jp>`#~f2Fbn#RhvbN~nA~_o0dZiirP8TA})q;J;E?f?|VvGD@g% zb9-Gow-;}my=daUBI3W2R;Yd?_^(tpplY`$-*NV3lu-G*^L{k(UlH+NNh{O=B>1mX zmY~?+fs7I=pUwx-#D7J^eIS_6aN(v|CO{t z9Y%uxN@X4OkAr?FqlC)G#4wupuZZ}sq!sEY68u*xOHiz8IHQCrP6?gPBWU8kBI3W2 zRwyIEf2Fbkm7fe%MhTUViBUB1UlH+NNh{RxB>1mXmY`VG35*h|_5qS7!-+KUUlH+N zNh{PzB>1mXf@+=L&iR3SGNXjb*Xk6S_^*igucQ^~bQ1hmDoaqTDl$r_{7GgkP5f6x z{8!Qnbv6n9E0qnXe67x5lu-Fvol6t{6%qfHv_f4#g8xco35r!+#3-Tir@o76;=dx| zzmis{%SiBFsjQ>^aXMeYD53I^dL>Q#S48|*(h4<>1pk%F5)`W%&nTht9|>Jc6aN(v z|CO{tO(4O4rLqLYv;77}2{mr=^g+vQ$akC@Y2v>k;=htss9Q+zU#VTNXfUlH+NNh{QyB>1mXmY~?+-HZ}yvEmPlFRa)-G^5>1#D68NQ1_GIzfxI(VuOVw66)ARXZ4;wm`{CEY2v>k;=hts zsFz6aU#SGuCcm9)8)Q|}7$sEg@0j4NeT63eDmLMEqCM3N?!a|CP!H zRNmkRj1sD5O2@ePAx->OMEqCM3iSyI{wtLwC|30uqlDU`wyKJKP80tX5&xC6LVZPo z|4JpO*2RVY-w#EcgWoVpsQd}}TblTAvINBje`S)Ow6Z=|B8tJN?M_QC&7QEvH_K!41Y38sQd}} zFPiwTS_{>jSfE;v;J;E?f?|WM7$sDGJGZ8Z|B8tJN?NU|9SQy`l_e-P*q%{Bow)e6 zzkYBUABq;HiT{d-|4Le+79+ubrLqB)pA1VdN~rt^c}bf1uZZ}sq!p?o3H~dUB`7wy z45Ng~Z|7xc;=dx|zmis{6-e-3sVqUU!A^`4ssmpJcV70fimgNw{}mDcm9#>2Cc%HD z5>%Ufi`HI2R<#5|B8tJN?M_IBEf&95|low)lm5r?8zvh@&47+{8uVVP^_vy zqlDTrJ&^aOiT{d-|4Le+4kW>Ur4p3R!P>Q|eM{xCAIvDBrgdI)I{#Gc5SsX}i1@Fh z6>1;}{wtLwC{{I?Q9}LpY{&KPy_Vld8A22P6%qfHv_cIf!GEQ)1l6``M=(mLdADzL z$DSRrCml%>{}mDcm9#<)C&7QEvH_KEhGQ5d)Zz0l8!`Vf9(|rB{wpH>D`|zQlHk8m zS%PAN$1+N&`NanBYqH~L;=dx|zmis{6G-r1sVqUU!7+>y%C~%}<>>Jo^noV+DsbZ(m#D7J^ewEj!{D8_u%z3@m~?~Ur8&}4J7!lRF|BI3W2R;Y<2 z_^(ug(v|wh+xd1z36r-n ze?`Q9C9P0XNbp~&EJ3lsrx_*GGI79ialo1QuV~`GBI3W2R;cGm@L#DcL9xLX7$wx< z<7~nKH{6H+iYERmBK|9Bg_=f!|4L;8D&JSHGD@iY9(;`^{wpH>D`|y#g9QJT$`TYC ze2Y;+RbKgN>sNl_o7dYk@m~?~Ur8&}yCnFpRFd?i%?6`Mxj)|Ex@m~?~ zUr8&}2PF8fRD#lHhQGg^S=C335^8obYAMuPuJWeJK^eZeT9_J3x@ z?h9_?Df%T%{8vQ$SJDdg4GI1$l?|xcWoK32F-oY*iofi@lVJ``{8vQ$SJDdgBMJU1 zl_e-vHJ4FB9Z*|U#eSxV|B8tJN?M_QCBc8C5>%VwYUxhTr}F|v36;;m-)Q2$BI3W2 zR;WKo@L#DcL9xMRZT=UP--FF*;=dx|zakc>RwVeZRFinhRWr`I<{kd2SX-L- zuZZ}sq}8h0liO;j!{AlUU1_f3vT2y!}2uoUlH+NNh?$*68u*x8&LV}yfUMN z%IDxJH1S^%@n1kS6{sBK|9Bh3Y|q z|4L;EidAjOD4{M*@4?Mz;=dx|zmis{ElKcSscb;yxAWGF5-Pt3x1ovuiirP8TA{Wh z!GEQ)1jVX$V3bh#d9@=={8vQ$SJDdAlLY^jN>J_cwW_TO)UJ#Y3a2DLSZ;dJ#D7J^ Qe?=TMfUh0=dtmqf0Z}4br2qf` literal 0 HcmV?d00001