From f77ea8a561d08e01c051a3d36a86d7fe356aa923 Mon Sep 17 00:00:00 2001 From: Lachlan Lindoy Date: Fri, 16 Aug 2024 11:03:35 +0100 Subject: [PATCH 1/2] Added HEOM code to repository with improved result IO and plotting scripts. --- engine/__pycache__/__init__.cpython-310.pyc | Bin 152 -> 165 bytes engine/heom/__init__.py | 3 + .../heom/__pycache__/__init__.cpython-310.pyc | Bin 0 -> 229 bytes engine/heom/__pycache__/aaa.cpython-310.pyc | Bin 0 -> 6213 bytes .../__pycache__/bath_data.cpython-310.pyc | Bin 0 -> 2773 bytes .../heom/__pycache__/heomif.cpython-310.pyc | Bin 0 -> 5779 bytes engine/heom/__pycache__/mps.cpython-310.pyc | Bin 0 -> 22713 bytes .../__pycache__/mps_utils.cpython-310.pyc | Bin 0 -> 5709 bytes engine/heom/aaa.py | 251 ++++++ engine/heom/bath_data.py | 102 +++ engine/heom/heomif.py | 175 ++++ engine/heom/mps.py | 780 ++++++++++++++++++ engine/heom/mps_utils.py | 212 +++++ heom_dynamics.py | 138 ++++ visualise_heom.py | 38 + 15 files changed, 1699 insertions(+) create mode 100644 engine/heom/__init__.py create mode 100644 engine/heom/__pycache__/__init__.cpython-310.pyc create mode 100644 engine/heom/__pycache__/aaa.cpython-310.pyc create mode 100644 engine/heom/__pycache__/bath_data.cpython-310.pyc create mode 100644 engine/heom/__pycache__/heomif.cpython-310.pyc create mode 100644 engine/heom/__pycache__/mps.cpython-310.pyc create mode 100644 engine/heom/__pycache__/mps_utils.cpython-310.pyc create mode 100644 engine/heom/aaa.py create mode 100644 engine/heom/bath_data.py create mode 100644 engine/heom/heomif.py create mode 100644 engine/heom/mps.py create mode 100644 engine/heom/mps_utils.py create mode 100644 heom_dynamics.py create mode 100644 visualise_heom.py diff --git a/engine/__pycache__/__init__.cpython-310.pyc b/engine/__pycache__/__init__.cpython-310.pyc index 0f30306ae33303f87cd8d6fa8dfc2409cace74b7..995913dd884e1a0245ea77b29cf6a50f7baa1331 100644 GIT binary patch delta 85 zcmbQixRjA6pO=@50SIPF?oXS@V{Pl9pOK%Ns-KgXoRO25r*C9nWD*~mmzkECniB6; nl$e`Zo?nz5?;jZNoS%{!pO}{t?~+)OsGpjbo|%_Av0M=VM=2a8 delta 72 zcmZ3=ID?TVpO=@50SF%a_D`M2V{KrqA6lGRRIHztnv+?PS(ccWl3JwilAm0fo0?Zr atnXNqoRL|Us-R)0sb5f;nx8wdMG*iu&=@5E diff --git a/engine/heom/__init__.py b/engine/heom/__init__.py new file mode 100644 index 0000000..2b7c0ec --- /dev/null +++ b/engine/heom/__init__.py @@ -0,0 +1,3 @@ +from .aaa import * +from .heomif import * +from .bath_data import * diff --git a/engine/heom/__pycache__/__init__.cpython-310.pyc b/engine/heom/__pycache__/__init__.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..a201734fb6c52a56c94be12537db4e4ef697ab11 GIT binary patch literal 229 zcmYk0K@I^y5JkIZ5*fk?xWHIw5laz4>=41aQgn~$nCX@=B97z~?$JB3as?|jSolfR z%TIoOm95n(jC8Nxq9O57n16X>1~kV8Ey*NL*y>CEqet|XEcBn8A|ZFYISYN~BcdzC z%BDJUopazO_fM+8z*Q)r`_qo+UQ6t36l&{b+*|a<_j2EsQ`_H mF;xg^P3&cvD!iO$e!!;y96TKm~FV5dx~Dw$|ITtCjZV zo{`qx>Pd>oU>V|60Yy>pBvCxze*mxi0-n2gfbftq@IV#t059e{J*#_MFjLdh_de(J zxqRKNzP_}DXRh_G|Ffqo>(5j;{Maa5!xwyvL|B5=trceeykW06jOuvZZ6sC_CT}C3 zT=7iasizw0m2@Msl0n%OiM!THRwRXof1gN+H2(b}BeM7ph*2>ha(CHEP7I16)C`JY zk;i|?9~L9)TrT)STNdVF%jPSLJpx76Qi172&@Rcqr(C_(EY+(8uAN&OC80e^i<^GA z>`>?4lfd7_H;*qkk3{j%+F^nzD{t+xycIc$DO)*X))rF^%C2&ktpZ>AH##4B{{AXf z;?1_U{oCzALc8IrAC@$en#(}D?N;3nwB6kdw6j^N2M;YdfYJZ?>g;N(;m_7f<<)wr zIeX>Ol^2TdG;5WbFN)V?so{58a=m!tt>S!3_{CCF6cFB8!6BAo4wHk(mT43B9%7u@sV-PVLag7jDqvsmyWY07 z0H$}Du&1oZR`z`!@_W`jc8^zV;RyGxClVqVT5EP5ViVp2MJZ8S(K2IopQoCQ&@c7` z7sQgAnX9+T*rBV)>{m6viX`5|>k1dGcA%j91p+97r-F|4t zbdeM!ZOtJHy}5d+aZ{A8j-dV@L|?zvz3{*3iw8*}N85GL2z}`dbQQ274 zUyxDu8Y~sc;f8FDD_8I_3prcZ4|tS-wIzfzM>S5|7N-@y-zOQe9aljaRZ>7VzQ^OV z!sjqI>Y%VxOJ$5nC|bwNzTcmjw&aedk}9#wBTsplGo`#;HULdS*;X8id~~jZ{wb9* z{ZqaEsW7c9~eLs2`Of3-FkTNZA}H9*a_;y8}rem#$h-245DvlaP|#Q+OKJC#y? z1?&CaL<6LJ_|lMEiku|I=Tv4V16k)LP-d3OqLkPxp>LMj1lzSE(1Q}%ObIli)YmKZ z9nd3ejL4_~a=cTP0u5wTY9KbNoj#Sts)JDVUW>5QfIw&@!<(?MtDjr126LD%8UmO7 zYXfSC!v6#Qo*fNG`Dg^R46lG!;#wZcq22bILMMb@_d{PkPZJXK&~~sXbgJBHu5DCdvN{!Fc+hSIzIK`|+0Zt5IdVOP zQ<`0uFVRTj-<&3<*AA8mBpJBo*WnXZ00gdI((G-`m#-V=VLXRs^Re?_HJVYN6`t}n z&EC;$S+j3xwxrnt0GyQCA-7??ruO#mfnhBc;WeVb1j)W5l6x8T!Pk+b*|bfd%VBAr zHh;JrPL8LaUYpW1pGMnDY!uEl$IigndT1ZF9mXH$Fd8o(h1-=rCOiS$_MI;2PtWOc z<8CoA? z&M`55*HbAmAtuGt1Lr*_%0$_y51!VAt@ep&pxH#8d|&^TjRVFh@-pzxM90Z_2^^$h z<>@VeSd}K|w`IdOb7qCqmO@*~w55)kxVji?b8J{fIibC~xGFpMYmqNDt(*1s-qxiCXvigHDptS#e%Ghc!pU z1yR6$j)KBHNc-tRz?x-a}MFY>Mwpy&u2=Y*BuPYK~zSzS!ABwE-(0 z8)8@BKf1Zl55fb-yx2Ly9$6g$Zo}RJGBJAw(BcrQ=)-GRK;E!*sDyfcxG%A0P}!r? zn9In!4`;4T^rvFA@p4b4Xo zDV^2#mAH$7VDOuCB>|PiwFucZD5jl+SEO1SwWnkH zB7P}=xM(H60Do(^#<6W2(h!c(4R`4?5P;{8uNvI)h*VILx)c&IOt8y1z}PsT(CtSD zL6wATpa6OpOCX`(3|iWdCCae?Sfz7ZO7>#s$~`m(bWWv!>LGM`Y?NJP-NEC#d*#T^ zl7^AI-lgGm=>RRafZ6RBvkTfwE^P2>PAoBRw_PJ#w9TzfH+;lV0xp7}aK@XHFRNa_5L~kR)lV4cWA!A#h(I!2sW5pr?oi zNcIj+GQ^0UwLY03d%1Jzr|4)XI(aUS8BxbGd z=+VvIyE3W$lYzvG^V5I*^u9%aG`Q}&{ zV{p0>(d_0U`pc_GYPsi*VJ@0y-^E@-K>=u#oSauyLJ$>(Vk%cBNix>kJTh8Th-2&tOg%uEcGkJc9|f-L8k) zY1dmwgVk0?{(x!|m0EQ}`Z`k*q8M!4q;{i$IvHRtOaz)6xWlmPt!jJ|^#;bp0dY*$ zDNo2!29$hA$w!oIQPM&JG=Ev!&6UeBHsvdu+9hm$S^ki^m<^jX{VJQ03MYv^l2q4Bm~pxOabJJY-Mg}}-7%ypl+ zdo-6r=Ij(73h69i(EC1JP24gUi-(lS2*}9Bn`^YhM%}*(6SIOdNIWYGp0dy^FX6df T+M7%#yca#!`p`?@ob&jZmC5~^V#$|$k%bRzLnCW%bI6dJsne%VM{#etP?=uoUrcJW zx>G!~KI$3N@=FNKwm7P&oQln{Yh3a}kd^ZY8{NIe_pI@XSL_h01iq_Pu6l7Yru{Ux zs~0_|63|(jV--8NbIQS;Q*eFJRsuH(T_f^=aeBqFandUT%P5dF4l@rlI8D_QFchPZxS|veJw5h2uw$zY=|%XPqpS(K!|O z)3rjaMi)Mc&J;3@;#@}W#yVa|^W`j07go}uzi>I$E0M(P{9xS#Wzo_7au6q}OEzmk z9Oj34OSpW7KQGh_CT@K9>oYf-^ZQRGqThktE>v?oE>?L72vcsk`ax^EqU*kVxytkWDd^c3OFh_%GHYyB)YV(FN1 zatOX;Z88@7=%cyiG*1B)98QOuyTJ3RJ8~;gxes9-u%gyMD>?#?a@*?@rQ`OmQas?N2?4+3kOD z#EqBE9RB0_iHA=mWZFp~W0L3WQ;(tpuQ8w-Pm%l^_YxjAs*Su+T?5gTw<#`(kw0p{ zcYf9AIM#2WFmB$DuNX%cy-+AZg*`V6j96S#2Wfl+!qnq9#&vt9)*0r>G4&{|JVq;g z1F&8-VqAJsd=Y)hH>HN_(-+iaa6>)G0|xNpCnLrk%i1FQJ>D-c#sK@1T~? zL8L@QY)Q!04h1=~E`1P>E2k_45fu^DOq+&1mtv%SP$VSzK3?#D#wr7P+ zQ(H`m;XpkL-G-kx?s>xhfryM+4hMrm=^czyXvA04gZNO$E(@{hhDQTC- z@BamSWA=`z=O#oM8FFijFi@}HQg<1Lc*PpW1X-EoVD&snjoZsgJr4GIE@P$QbyKrV=%_dD)j(Xe z%V+9M49WH-X%1k;0QBVvk~3mfv_NQ{aBpoK-nHcY-aB^>)gg87#{W+kRBVe~1C=4( zVLAeD*BZ^T0a%>)0CUKK1mmtJ#SiWXQ0~?x8*+7odX*=gx5;F_om-$D=3LBi^(jV; z5n$hKSjW62YhagPt0}#zKF4I;O8=TO@)0!tCI{P4#`cjZ0tAXD8_*M=4bV2wHqm-; zOsfs@A)c|AM}p;xhD-H4&FlhNlpX88?9bb`Z>w1>W?Nrz$|#+2WO<8gQ76-3%i7!0 z$v7x1#0kaYK)IuG40{_CHf#AkIGa*v*t8Bc)5TB^hMKhMEMK;SXy_PtN4-D_UxYA? zaVpLdZ9#<0_PKEu*J2Wv_l8vVMJ>>sC>vPjX@~qA{+1w+bNr)-A`?VWzmP)^NMu@3^w}`( zjaSrLw1;4-&Xd^1ukBG%R`U?VF*XPLfU!R{eZbDbd&Hj?5m5WlVn`gTo~C}o)be3} wus#kPhYUx2F#fMV;b0Ta*jaSK{~c}ox5w8$Nwi}B2>7f;HVXWD;LisC0cjgt$N&HU literal 0 HcmV?d00001 diff --git a/engine/heom/__pycache__/heomif.cpython-310.pyc b/engine/heom/__pycache__/heomif.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..01bab88489fbd1dc27c6ac4f7db15eee9c36765f GIT binary patch literal 5779 zcmb_gU2hy$8J=@yW_Nb|5hvNij#HBKgSJcCI0=Q864IneLx99>3}+j(x^JF#?=H$d(3jzbcBh-%MV~d@RXxdcJWj3R_{I8W*+S32w%E zvGLj)jq|Yz8h)f27yQ(p3Zl6%3Z~|RxMM`vI^JE$tIfF6Sx)^l)*dGRTtb`@dJMP! z=Vr&TS?yv{U7{GK;ART32c>+Ek4-nWUuTHjdnxuoe2PD=zrNfgURxL8@2fzel3EK; zL1wi8`+`(rP4=AB&R~lZ)~dK*{cybJLK#lxDtjF2VHtln6Ht|?Ss8lJC&w|vX=`aL zEoE|5^hy|0&Po~721{|Kp$}K3E8RfT2sYpCvzP0h=G)}%?Lfzg!M58DZk;^!V(w}` znhSC#S<-24`?HCG*kE2S=Wa9`CqXXNVlHNKk>qmrWG=$oogp5(k+Lr@b_rQR*Q|*# z`K<8du+WdAdbLTEx2_5beSH3a$^-DhZhNK7PR5jteGX!ExQ9I(w93jMTEGYYotPcA zfLNJ`Wo3ieL1p#4%+}{I4(J)=qb%DqarXh>ceqZ%_oxX+I>HfS;(1e=6q8iMQ|M}v8MiK$`67lYo> z@-}I~SfP!&Lo3W=&Co`zpof2LPVshe|Ll{xTkN{Ox!h5qKbMz$V_|hWmlsZ%P3H3U z;(c!)#oaolO6ccL3Vj^69Tr1t;ne0$e~ALL=_+_Y=IAF5|8WYrk1 zPtN2+2}jTh9Py%X?C;)lgl0Hsa-pq@7D(zdxc(Le;7UfZ79D*S*pz@I(94`v7wz_f zVE8b1dC024qf)tG-$4Y;TzFxbHUca?Jcs|f)2n0^L|p4Vsem)(cDV(8n6a1#qao>C z)W_hO_s?&b&Ul*fKux;=q}-YTGM7m=FSWx6u`+jpTis4x)`1CGaLsKu+YPKr5$sy} zB*e~qPGFv}~8FQ7SLAaaBU537k|;snU9CiFMZxWj7cTYam|Q6bO4$7Ci}!B}Lc z)P{{p*dl^2E!5l%w%~i`G@iJ_$Nxr+iA z#4w?ZW_N37A3{@iISm>r>;zE~#!>PLCP=9GNihLUQETMJeZnb^{pVOwlKKJ$?C=c| zGcy9?8y}<4w@t=ru?CNrqX%s2(`e%Buy{nBtO6Z?V-8upYUbUsQzxs!sw>bh`lw~) zWVUh9V@qr0b!4Xx84K;+KXz$v<^gwa8f0@8Z?_Hq)j3pookiD+xL*W2b9+T3 z`X)-?2rL_zP6Yc8^d+=6$%$JRZKfZQzK;s{NP_2WU~1j5)~y>@m3b48 zm0+E%wh6S^XtT-wgpqg;Yy8?2g4ro6Q>P^(NAhI$rD=!+C4>I#*M+h*gS?tWig zA_Wds67#xsCwFH3bY2^b`eSxuQM>f&YgZZ^IDp4nSg23S!=MlT5$>*_{XQ6&Fd5UM z`O@@Ku`2Zem>i%$%2?^yNMlIOXh%4cOOk?Oj5HpV|6n{x62n!imiPc%j(M%yEt%<+1j?2LyR$mSVhuNI9F9xg>z5Opa)5$ zAKnO>j51pnIkB&0p3ty8Fd{bf#GAW8rbTrxkW4T7R38?obysQMPA3b_y2UbMq-WZ53)wlh!|gxclFG_wGRjK`TX~T?wM)4nudKjdza1JKHOb z^k&>h!ZgSyi~4MgQ&S37gKO~smE82Z`Q}&uu|J8P6qT^)a~bBjTVOmJ^#LM{0$R9rjBF_tT0$qf$6(kv{Zj zwNM0u6dcxG02QV<_?SZx5)PxHJun0p;jXlqa>z+@2r2^G(ga-DQ3BzyjOUSr@Tp=B zwBu;CP*Luh14A+d0ZK;NqQqwUAzC(`tSY`o873iN5+2gY63hssp~?uAuV>{uL+chO znRMY5=i!0iHjYifVeIC7sK0}9?u1d8>H{d6qd?s*QWO$eaV!|l%dNJbLVZQPM@0{u z6i&1HCWTG$5=H;?kdMR~mNge6sR@$+N_F&nd?I(Kp}Z8e!+7~dUYceSn3(q|Sh~qs zW+`_2a1x};-G=dGU&YB!F(%=1qu``b5cC8*`ItQ^_o3#Isypc1ar7oX^>zF^R6x}@ zI;2H{gY;n)M^qsDW#}6ml{yRVs!WbD{RYPPAsxi8McNiiL%Z5dzv*LgBE~kykodQf z$Rd@bG~*nsK<5<7o`Y@8iIxNz*gcR4)+y%0)6gbT6`qd^1u<#C#}voE6(mT!-(m+3 zQpOSWYb-Hg&dZVVwf0x?O5uPiP7zfNw&s(7$w$PMTwV!xgFO`7O10li8_?NrFjvAa z2aY%hs3&kTtBNBMd5O>onm>K%=O_O5(aNi5?*0RI@af)vKKIVAf3@_j9rQ`m>0vaK z&)ggXp`7!x`kh(13C!L#LupaQrmK1M0=Ay4ki`*B|v_!9Fp2ElRokUj&E z+eqK}k?H0JLU3z2YUUy>(B+4Q+&kU&JG06^^Ly0HH=<+K?b>;L{B%2R`t9Tl3dI0T z0^{eAcG!WN(knNl{s|{3bz@|t$wopC9d|&7_lpm^5eG|Q3gr&a4c?Ylp_jLR@ zy~5_(6lP)L@myX17~+pCybS?-!iV6a-{By`SgS0(NAWWXGuPmDLVmT;A5r`Q0(mC2 zT$fkPT~C|JR$SCqQJqHmp;Hu$Dhda@!C~Y0?&SMgvWpgtz`5M%t`wNNNi#&5nkl|a z(RFqopUu&aiENSNUSf$UmB~L$a5C&={L5f(-g&@F`%uB6K%vX+;0%p7A;!ZNk~Ewd fJlnIYPPOWdd63ni>V9v?t5k=*nzyeyI%@w9FQkz@ literal 0 HcmV?d00001 diff --git a/engine/heom/__pycache__/mps.cpython-310.pyc b/engine/heom/__pycache__/mps.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..1715a82e05b9bf41c6f752014c1d3931f3fc0116 GIT binary patch literal 22713 zcmd^nd5|2}d0$^M)6+Y9VzF3&AUFgF5KDp!5~N9xf=KX!M2cJv2#Rc%WV=dRZ@W(dfR>@uk7>I z%<#MD!J2#5Mml3B_B{)H?)~%ea28iI;#gM5ajZ8SWy`;#oX}nIl&ic&U**D_%7?i{ zPZgB^b`B%-mgiJa1sJ!`Y2m6NRl*xz4XY8%t*A!T818`@R};7osXb~E_mbMHrf?rt z`_we;y;XbbJP-WaF)Sc=u?t9b`brZm$~`(AaoI*$94 zno;-QzE9n&PT)SR?o+e4?^h?){kR`c52#bP-=-c^58-}LZK&^24`ZgcBln0pjocyi zsCo=#cc{zijQTL%l-1+v3B0)zbRIIOQdiVD^%3Nbs;lZb z^*nON)W_5dXnD7KQO)6gT)m`T#(hS;qF%-Q9(7(_!2Mo~^>OtYTAV=J*VXqTcb~eZ zKB3+~ZdQF#y@^^U@#Zb{DZIH~HPxq86}bo0XVijH$eltiREx+xh}@D|M(!cxYH9_! z?@<@kD*AX>t*JWhkEn)P$NjXrq%`i2;yWT0qvm4*J*~uDTaV^6Aj3*>%}_U@3yHU| zTw9e?FBDOb{S6_*=cjQ+F9YCO`Ic3&V!LI>&Wam*i_UdNm+|cMKDVruTnq17UMq(d zg_f`E8wJ~{Uze77w0Ulz&2w>~R`49rDK zzZ^lMNBcTZY@e0&n)c|R#~__!ux*YHV07!EqpvFq-&#~2(m?VoT3-kgr@o#P>Z+=B^;+Vs)}lBWpWmphs>-5ntX1TFl(@BeoVcrDJ@Fd#FiQMN z9M+?Tj*?=f7F8}+S8FO6t;{zz*wwYlU3*Vl7mabGo_~UTNrym6u+rJl#-XrCL{&XR2}aWLRIS)x(p^VPowiKzerl z8a}APm8(>^P|Mz`t)#$-QJx@X5*B{Ks2o!k)GRr&(ebS@Q+9&;E!6+#x}N zvp@o6wS8R0*a8{|A`lItvDi}pxX;^o@ASMsi1z?(ya!SM#C_Ze>P4W1yU%JBFFMgI zz{X=8Wi8pQ;v_+@&?+kTqNke%uv(7F?Z=2r7u^=;ERa*l9)*AKl<*Ke#=U2mWEh<;*4KP$*l|DiDq|ZqDGabjI z#N~y?+IlsvVOq1A+*T75-hCf`7iQc6;09KfSY?#SFq07^N$~={;90F3T2Jt~^y)@E z#>~yDN!H8peNI>FOQGJ!7kimZ@$Ha6GnQf_jAn-Peij^HavPHl*%+8g;>{5mojG7* zO;F4vWi}c^awo3n1QO5o>@i>gFW>R)QhFT*mYlHf!rMc3UT|g7Ztf+9Y$GZQ9{S1t zh&H?gRsJEC%%B3O@mBy^HmK@;00coaXnQflkd{YIp`05|+q+aEaG;)xdOr|@)2skZ zQ1@ip1ERSr1-yT@#diJ*=(5Vi?j~sKN^Z+m`9*;44R>>h^|tKI(w4P3Y|0+7(D&iq zc6`s;X^E@99OvVambE#$WnBk~V$9?19L665nn4J=;n)_&8pk!kS+#(I9AV4udve+` z=T0M#)&ke34ZZl@NB#>HIp+pU7~E?WBA@oj2I$W1A$ zgZZIs(ZyVUw&m)7gE0?oSq5RV8P5#Q3A73xBzb_WYOEQ|pg8d^6G)_9RBM8Vt%i&7 zO*?T_P2ZY$E=6bpH1G>->ehqL)m3=4+}McAjm5Gsov6In&}HznI4qO-;DHTlU`!DZ zu~A>b2>AQMCuVLJbW7YQ)>_a=aL(XTv|I%f^>uhT)KRFCfQY;rUVY%yX?+ixC;56q zuT>%Otm}qKyvw!lN>W%{t;XPoNkNCwfj(-gl6FK%5G~ghWAMhh!Y4IiNK0UX_|{80 zHfT$HH|Q^zo(Exz!DNHih`t8j$Jt{b-z_haVZ+}Wb$($u;~Nww3iYr}`$;{+h8|~4 zCawmN`WT8oj4OH^NdUABfWY4SzB_{4Y5RcVIg|E?eGIHGaJ=`t!k8~;?;NyC_CeI2 zay;a`Kk|l}qeSx_y0gGt&Hk7t*>Q12D@Xvu;DJE!yzIU5-f!=Hgc^}nLQ&sBYGT}tp2UN8Dv-W(gwI$EIezbWWX~Qt=1Mo14#vJ z93CQJ1X@5LzgAs`Os1zXpkW%^Km*T5kh=+j$rs#yv**tFEf$F>aJKxx|ENrR2r!8 z9@i|en`50dOA8k`K=#M|B8tj@KZ=au51{i-6vVC(pIQ#oQwKuPp6rbP8A0bBg3e|k zhJcZQ)?0KywNYR7Ac9IV?xqk_Ll&l zUt-c>`iY-O;WJ)$;}iR`evb79_5c77WuO_MJUWDg99#&1GX;cp-VdA++#&CRiEt-gSyIb;w(-tImYDLlvOj0B=2ps~Xs`YbAc6_?N@-);sR z!f=<}N8r+#5tKcRrJHe#&aW#d&^rzDAr$;D3*SPL4kEm7+aP^F62lw-ijPm>f3|S0 z;e-PCyM0(P6kQ3!b!?al_18J_Cpoe(pq`Sx;RzZY@OTIhv31=7Ze!P4c7=B^p9<-L zi@S?EWXB5SrCbT5tpaH?n_^z`I}9=kKvCp2!Z+#hjG?`WtHXEzwl%5rkKu~Q^QlIf z{{Z52)93C`Ht3r-!J3zVH906|)C17M0Cw1vq*RN#%Jw^GCqR%-v8cCj8)fhC>DSOq zh~l4)Pv5zlA;WIgyVEsn@ph0cP#xWToExSNckyA zXtUUv#qMkU1{(bt)_P{QwN~zafC<)HtOGC?^+6ReAzN>~fCm9o2;HM1sZ=`p^{%V@ zi)^HjH1Fou4K&+vwK=FPxel4zX$>i+{{pV}de;yVr_>S7qh|DRZu#`M3EWYR6G);f`6)HmVmiG3R3suosg#8;f@Z1fO29(iu zVN|o&S9yW$f+3iY!HBrsB=Lcfw7XHe8O`#vo9^gvv+~ofZ_yWfZW@{k4 z`?4`069T@$p9XxPL5`5`d;sYE7ME3^_ZYvnztxVwdoNqOfCT0g0_TqMS$E&a6MOb; zSE8|wdcVQC%Sf6N+b7e9wk|615*Q~hQFfbFLu=Dc)%?jeS(wvu0NpT#({=~=KG^@R zZC}DNgSESJ20Hn10H&k3>He9g(puuos~%-S0IZDK$<$i)YHe*}P4wID!v6!52_2-F zXcCw6IdC|EteqR|j>^A{t~#hZKY+^CCQuMo9H2240eM>O@eJsJ9#D>R27ZX%-I`dC zPxrs+L3e>amBYUvJkTB91lW9t?kIz$#E=E#Io5tFB@AdWSeNBJgj5oWLCo3mjwa63zB#LFXsB11Tx24UnR2oLC={-B z7b}r^=rtxnQ3#ZMBt$HPb(RlMkfD?uaAw~|xv0$~8)XPE6Zk)FKV*CV!_OHRtfzJ# z_>w#qa0&QwU$HCJmIW^p$+MG!@C8a`n~=iHz7nAM)ZJOp+hFU9J{5b|M;AmL_P~~d zeSL~|>}D>`t*6+7j@u)?>Eab))>K zO0!oK{&TdA4kGcLG3?U>_KDh^Z{GuEdOy$+BIh3&1t)9Lw^gJv6cA*qsHI`SM4O7B z7E=Hf$|Lgo3M4u)F$(Pgnt@ABr$~sc!!tJ06NEgfa-coqd>Rw<(=q}Ebt$oUc_br%o!g1G-Zg_&gba2>8Nz#CFuL&OeB&o zo<(cwCmC4E_Z3B5ASePDY1+L6KMcOaBM3n8*Y}~7W>cZ6_zLP7J=;wWzERq{vQLMU zY7Tr0S41=+&4g8CoD8AQ)YnHS`*w+1Jm!m?*|Kf3UzB?wxpe#|nB!s8lIg*&<;IQG zBv&U^gLm#RhRUAN;%hfkW&IjkHJJ=h20zvt^#zz+-$!9Y87u(3H_MzAE%0^rN1*Be z)>rXB+Ky?Hfl)XE^Aa!W7~Vl(!-vQ_*q)Ea9vxwKQTAt_WItm zzNbHo+M?BoYs=Jyueek9!L+8U&k?6dn`I-x6e&J7&yfQ6tw!0nGRw==NZdU0VOTHA zs)ec?H-uy_z+nReba}l2D=rMrvol3wtT~H5&E|8XyP%O7h-vkPVk^dCU0YhlgK_zY z!56^PHRyGBV}o7wZ;Ad4-nqb|h*W}hT%bOzxZuI>2ksc*3Z@?K{pLQx)((5#hZTgd zj#S;k1N@iV0@4Ck-^*OviS12Zzy(OekU2JUsf(lgWVRO1;^MG3VUEjO9F#=hatFkw z003|3&b4mB6hQ(3bGPW-!p~^F7QNaqwMY5*1QL<9V_jW{GmF!+eCsq81{N$S8lyg_ z<|LPml*7ZaA^MiM9YRI0!IthFR${VmmtOZ9!Jdr0h1-01;teZn77>MVKF_zn+(SnhFaKprZN?Jnk|=~N&^(c)GjCtTx~4XvGoW@7&OM-g&Aqu zK#aYCGE~h-5r79xkPhnr*sSKAj%{uZ-?s zFEEVSVd%jK*cUu9u*V<5(-`gYV#70ywrOKTRsliMFS02B!EM-lBC8Yal^mrc(QMCS*~6V89T8fuS09$5kh!x-TRV+6$v!Q*UC>zXis^;s4w$S zSjfJBcw}Z~?+FT!B2MC_JA4bsB9gFRn2>T2jmTdHH|_iI?m2GS2r?94g`eSm&~A=& zIT#20-?3*MNI!^uGenw9jr#`e5_<;2Pr7IMIemi@{tA=7z(nAjG7Y%(!1sw(z3WGQ z7HwV~{Kye}kzq}tbAJUHeyK=|#U5<5=m*pjSPkllcSs#?pJm()C^>?X9Pel&MmvZ# zj3?m=oiw$n6V{2bEXg541GWT21Quc=<`}k$*awryM=(*8~&SVoihmdLJ4Y} zv05;HLq}|0APPnc%WE|SL%S6TUB1CB}W?6@qHn6h5UDPEL1Qr6jYITr-wc1r6 z7*MWJuQ9lZIG;ut3E3610wJ;a{B|nQH7N9td<6smu}f8Ob&QyXR3Nw_9e)-W+(F(X zk07D7S{*1@B?-y z8}J|&3BUuX1EN($%*qKoI06s(mdhAXs9*k@G|Hh4pn!Kk0|Z{iw~F?19X(X?qI*C= z6PhBM_(PZnMBOeWc?mCs)g>c_Am9td@Mc6`?nbh2V_iWVtFTB2dK!P3g@}BOeGK3v zx54<)euTLJwn?&brI$F3IG*d<55rO1OM$6H_DLkf+K#p7H}Q}O#W4)&@Uqi-1MZ2H zhPnWObP}PA!cBpT4iy5v7sdMt7~^4YYr(3A_$1nT#A^>>HlT?U@kytEcP@0qhY|E4 zvv9t!Hn-W6o8B7AIe?J#Vb_lCt z=?9QB@0g?guHM*KT9)WPgvr+9o~Wdm;X&zglx>MN`Ye+tnH*&z<~uRa(cqz6Z=JD{y4|d{tW|gFb)a`C@y?D#E2)X+d8>m@5!Lm!aU{vvK7 zU1vHPKTdyx$-qbDWOiTUw+OO={}bm^-x@!-PjBT&oe6XGu*qZ-t{zs@rbLHbbf##@ zA^_OT!5-c`408Kes%u1#m5)Z{Km?;y_s1L692=|ySC5cF*|q0)m_6M6h&QWLzQOsE z@C2fVXHI`9>sxRxl7~wc{4_G1`HyF_r;gmRAH`>qtWXc+Mhw2ye-4fHk0Ze(NrMyn z>fqWk`gv_*wNiNp(|eC}kVYKOop>K)G>dPR3V*gano!2UIdA7i~1xZE;27!z@I~g3#cKnZ^9@BxJ6tjF@Mt&yA$3SXJBSg zqM{7GPx?(4p{wtT3mwVc7eKln(T?|{2b?aT>I`ZV;6mgOMuJTWIAAYfLzZkh(W_oXII!D8yTtKR8Vn;7K64f`=qkQKfeet-zLXS2-+<&2A}ePB z5Zag;s}26nu)WlE)QgFSnM8vdI|9{t^;!jASmC$)61G+dRbDi+z$sBXvqzwxopJ)C zFj@(r>Tk9{CX@Y zH7PDYg1mSI<``I~a}wTTz zwBof5@J0!OqPr3&SWF&LxM9Dfx=9L|FTVUzrH->BM!p(mKR?UjVL0!_v}nS5Co;)2 zti)Fu0?RWKM&K1zCRW;hnG0XHPV2j~R?~l(T{pwF2a zmkk05gDLaXf!){lw}?z9P>D7qpTT@i0I4m2J7fY@U`olxH4@mt z#r+RlUYaCGou1`t#1DCQ46m<|FQ6NS(MIe$Omq;ZD~jug(pNox@V4M(GJ!U9lH z1o^{Z)*B?LAgJvUNam1mV-Q6UQ=F@#EDydp`~_zbed6Gfx9NlE${U7)w})Cd7sB#X z-e;{6NFq(+UqEjqNcSVHpnq1a;YC+YXB0(eA3^DG?>sxDzVPp$#j_>XThQ`(*3Se2 zVT=H8zlrmELoxG|@3F-rf-F?{s)tJD^|V?fWrIgHd21u{;8-9r~i63<8TFR>`Cap>4n_01{|Cg0N$d>sOb$%Bd(h2znM%z;2XL$YEC-1u0xh>{k-|eB!*e zi~uK$rhk-^7=U0R_boy&)e{D?lf)f!ra}59LG*y3cs>L}x{)M&5W_?GH%}$EIYrX9 zosh#lzuxGxEE0D+bnkO`fKUNhQdsT281yfO_|O6^RDBq$a$8Q@t&jE9$oJy3snK#* zTu7lV1P<7efFICo1zQl?E>E;H;E1>23}6+M01+ks2=nL(D}t4pq7Ak9l4~QH(PiHZ|c6OPv$J@bl(VsCVXWdP#mQ#5q^hwNB~3 zfiB?9E>#R~6z-lGKLqAPF#UDZ(!a~(8WVACv*o~7Wr9G)*6-m564ws(eO6Z7CyG)T z?1-X&AB){R=+bgTbS!{o!o3$m=)2A^t-YDeIeQPQ6UtHMilx%{5o`o@(OJklU;$+h zpq;o0Rle=ZDDY=@1XR#e_JXr%LvK2yrY#==WUe)|#S_y0rL&-3X`krpb@WxFr2_(G zURUA*sH*5oct)%XDA6)n2lUMIlmw_!Q3T&o8$xbef)G2b0#Lgl2)-l03d|H|rX5wn zeH7szW1x=1=;dEvuH!K%IYirhDhNN(o@h;gi$MHEjbDjJRz_PRe8SjctqBOlff^B^ z7*C+TQifH4*xVj(jj^Bh*rl&#wSh-_*hf4L!FUhXs<(ar6C#1o4yW`+smq~<33K{d zjb*_TSeJ3J08MNUr6nVzs4Uuxgf1gi1e=HxPTbn-PH*_KV@N7B6<9z%ikp6iN!M;r z=HthZ%x2(rz-rz_`qT^5J4w*ZspsTuoE(&c)woe!1u4f#*gJdQldD%js&S5XV;$#= zU_6OZ6Eb0&mFeI$J|rHcGagrC^5+@7*$=8L=%44O?C@oeVj!a@oEaF5+AZZ^bJ$PJ z(tYk;*py+G#+0A}&+5O*{(q6lVJ5%BklGm||W@o%&8PK4=iq9F0sptBebDP-&ySepb|{~acO zmx++J-#{*L>xkMs$9G+V-BBL&KS9BL1h8k2nV@%B#MsGHGoZC_1o=D!Txt%~8puMa zO?YtAln_?M=#r9i0!(z$_1)AW(?fUr=pkughrTc30SMyA5kQ`xg2dQ8FC-<1>(eJ0 z;fLR1pu&DQt3^JToh*V5(a|8FOYoN;QvQ$~jPJtVR|OkfkD)-GH6b}XU9uaSLFtjLy;&P|lB z(j-rD97|63*du-Q@t#C9YtvWP)5{2lWFH?v4-i}0E>83F`6WC^ont8cn&c<>&@fY9^z;ueP)ahA!(4s;#s?ldI4vk z!6!MUQIDM5<@2;S7MVSpv z{QV8RPCz<_wsec(-IU zoG8h7L$izB{6#l9B4tqOM9E77Oqhp|aW`U&(x8ZL8wvCUyqXqp`5+FAA$){I{yqgD zNPgc4)-@a$Xn9c6aHNxWTG3#ML)}fk+v++M$~e2o|U=d z3DBJ?m7-eOEg!0-3o*ld-L@JKo5VP5YQ+J5R-k9L5V!8*&(C0nMKewdP)~6Q<4E7> zXg%X}Yp`!VWy)*!Sq|ENP3F(<z&Jh7uY5~7|4t2=GoG@3Q`MMXS%4xDq%Ioa zrR1Qup6dFLt4ayiG~_ZWU#Z2*_#uwP#ZUts5WEiNSOebB+SmPo4hR+iNvls_aPfki z@jR~g$Vx?If|s;Se-ik>zU`;5-L9mR*w+TcbP<-migDiKSkrLFiLg|`ly^RTeePy{ z_us=H%-2VoGshW@>~Os zblnZ$PV)0|zEyrXDKfmf%8JfHv{P7f0{@YG%NOx|^7Aw!@Nv`RgwSrM3~n(61NmVD zEaDEuiyWCx@Q`>DqW64 zA*1K9%tz`eKaXDYzeOfPyl7*77sS*5j*b33lmE!1+j$v}1axWplYcUxAGfF6odI21 z_W?TPPF`YH_!ct#EKH&eXt{)=N3`UpVapmh>gUM!Kc?%*2q z?jfdNv{R4K2Dt$*VBUBE!J60;*2PGq6K?)+^AV!4T-a z=^@nIJ!1<0W#&k-GqaBr(Yj+)A+Z8UA{L+2i6Sy_Cqb0XyKF$(Z1k72%njIv9^pGu zL5RlB;Fj1=@vO5UT#54}en9A{TH@iS4%WBrZUaB8vRYmL7I&9~AI=G?75sOOJICBZ zqK0H-E3MsUee2uJ-$xg+o6sYSK+@S;hvEq>bBxy#FRu-j{5+#aWG|oP2_ak%1zHen zvC|Z9`}UNlhUlm0)WpwkNiFVZXXkl747J=#Q*mz7hY(EYL0?-gTjC8_x%7@V7!hxZ zB9NFmrX7BRVkJm-dX4sD*=bt&hkBzXewwSEn)SHXJ+{5hPD6l{4c%e?-=*oYIfWm= znpcY_%JVC0<@+8l<8f_8EjFjq%!#sExmeFi>Wdep#Tqgv%8M5poq|R-QU`v_z@0}# zRCmEi43wf1pXJ2+VJgXoh^yXM(}S>OyfOcUE#6}yO79;dm*favrv?FRL@eD`E1{9m z{~6(op)X{V0%^`Y_Hk#zr9GGmFCa}k%Kgpy0FKB&O-snf*O7rdXl+>fF(4%G_W-=C zE!ZEh?TEktx!_L`;(`H}t^^leTDZJS3iNj9ya5ZyQzF2W zKe%SXr|v|1`NcAQjJa2soM-Y0CT}p2MfntSXPJD6$)}mf%FHuIB4;GJ26O97 zG?R$QWhSpPA!z6}lg}`Do5`PN@?|F9$K6wk-(>P@OvIA(Pno0fM)XDfFPZy8CjXg<>?Ze6_NQmhgwd8+kU?~`iYMPX=1lFM zIOreshkehVgkPllUqrr`{ufX_;SUA5!1Kq0l0O}c2gM-3HGz_Gzu@QnLf{9J!BM{y zc)^iiIBQ!Z1c-!VF4?bbIEEeZyI?NufMw5@ zl~E+$?YYRjl|A)gF-p3}&S+@K27=KEIA(sJCLKyJz%@gCgS`p;=@(QqD;3bIcmqak zdER}*F&U9jXT(@IF-Y2-Ied`x{ohfHJU2$R4H(h=P!;MaaKYofo4+)M$E zm5sIaYX&9wyAzd-xV9Q~*H9+9D{LhG8G~bwladIF_7Rjl;#r^V{sT&}by>uGeC}Fh&7m5P@BAAST3d_)e%`CD_Ix6hs-TN2}SDX5^k; zd%e}VWKF`wsp6KZTx=>A7l-_Vmh~@c9RAs8yoM*+MIkJ~`qoEGaADuFKH}npaD{h^ePoN0@X_Ol z$3!41=y63=)X?LJIWdoS2_p;Q7)E??QY?z&sFg)sJc6-6EQw{jE8=mnB33a{MQu$y zids#uOIEXSAE#Tb4yV`nM|Gm^O+48PC^V?swxq8u&CXh*aXVX}_O>gRlqI;fmaKjD zd!{YyYPTfAy1mCQ_o7Ceh-kO*%u}M-h_gnzqZ(xaq4nHSEy)czoeyKyAi@5pcv-ALtj`iwRcq9&(?=Gb8iq1<3!QbGLntYe73-<%we*G=H(8ZI{8r# zJ%$Mq6>V*S4Xf73)s8%?xv;i*c1*jPZ#(js=ICdp?|H@UfE}9QZENc;t6SSBN1pb; zMlPS$F6yrKuEv$K3x`q8lVR?Ka&r*w=1z#R zS<2a^+*WBncay7WB66?OiwAL%`*GGzrRt@*lVsP?pD4MLbPQ*U!o#KVx`JD4_4*SX6! z!NxDxT@F3rt)^o%H>}k>Nqp&e?^Wix=^nnZ3=b zx{N2Pbh$7*co|>-e zS{Izf38mTj;CqSqK3#|1_|bJuzxB@Ldpxgocl!N;iJd6RgHAda?tra3*Y7b&oRTD6 zxrE|gAwg4@%hdY_il%SaEaCaALY7`XPQw0XZbKW5=H!lu!%ewDv*+X*O?;G!6DV>w zl_HYnsNB9VP&vOa$l3d{MuYrHPcBl0UCjAKmD^W(YE~7hQ&kwIPh+xC=qg)cE4<3; ze38|m%syMh-*-;_7bTM$n0G{>%jh6GA*IazjDkdci&TT{v$x>zMz&%QkS&~BmT)0( zTX-}>+h){(`<}MkKU!TTeDwI*HRB%5`=mB1X%|LP66IA(*`HzGWvB(Su05-qEf)sl zZCamMUod(uVnn+lz*#=2nfAIiOsUf9koio0M6iSB5}u4$ya{do#Q&I&*a-Tf*(U#- zUxTJh`3B44>n})JnX}z#zch{_IMRN&W+cY!L_EmW&-K$z z*w4H<^5KQZFbmMEK|W;ozzW&%c2|LxOMZoq9M4 z#R4FWy&c^=LI=?>;n7#*%M6da5iCqtt5(4uP(t(?BJOcXaG2Y;&2Dc~>lG9a-0L_7vcp7w2A!M`JteMLkZiblkHKe- z2`6@%Q+_r%r7&g&T-l8>;7U14b_S6QRg`;hpeh1}gu`J3jIsMsVvO1tud(rF*;p`G zlL+}L_RL+>23hU^az*ZvjHrPwdDaYHAjez~`Xi#j2zc?^vivayGqUzYw&2aNIljmN zM(iXHm?Ywlgf67K2Tq*ybYicHDJcCu{*#)pi=Mi*&-b~pMy|-6!^M;9hE2H-u_^Lb zMqkRLFEC<}K?N&DwjqFr@;w2*H3<53(ua zXC<+x>E`xcJj`7aRms=TH5CuAlxWSF$ocN9+UCa}aM-UfJtM@+0?O z!9&`f{Q<^N0OdX(Ih%-%1OX$sjo63X1-}J94#-2VuN(j)^{^FqII!%mqcS`g)9`vU zt8O7KDoxvK1}e9_t&(=`_!?sbyX0%(X$fqDQDszBWla$2!5`LiZ~(ToG1vy6B+wPW zD|!srR&=#!9Rk}b_NW4E9bE$rT$BK~BlzTs4si0A#)2{cEWnq%bFxG;uz&TtW&v`2 zf?S#ni08{z!2#_4_wgxzGpEq~YH&N?L{SQM7pLHi`Pp$n?g{gXanS2uJ3-N6GiVh~ zh2Zy{%kpu&8Y*pcQfN`wQ4NG!jU-K; zN5<7iuWm&hwZ8Y02h5v}H*UnL*AQ{H8%Z#$0l2vuNlFeWl*rbL6!OAD_sC{9F`-AJ z7x%H1YGlK(6XA#fpaywOw%+veAkN|>Q(@AH^59C?--%4}i3G_Uo6nmq+6Rk}d-dZ^ zWK08@7<)ZuZ{+Ul>L&6r_YHc@D)B<67j?GV6e}T}0;BdJpN!llg`5j`o`JqP-^kyg z>swN#_`yi^B~;+pNK>nb*dRF>Xyoh~(1d=D_6kLEcrQ`bZTbne>=uw9-t4K|=||ln zG6v;5?&Fa;Pdm9SVmV9p^3G5I?V*0-%(2~pSY-4k&|v_$$OOseO#VlvJF7)Lg&b%E2ab%#rG%K0^?JDNKL4ReQ}p=7`SJq8+( zJpmffIs^?Sr!sk1aFljk!3P1X>)NoJie?msiF*S54*dtm#zQH{R0Dy75K%9uHi;>Z^*8Y z;8=d(0o`AfMRgyrWjsZ#B}1q}5kH%Fm!OF&EKqTb zimNETuD_cgaEb*?SvK(f1A%Fxnt`PBY?@*k;L7te!AJd4ztmmjU16-e=FJ~X*?7dk#suvEV6T%(FNQBzmR1W@v=^4ev zK0xXuKBp{DCTcGM)QU{>6!-1vhzAo}vDf$n=uYvfEhUy8F*b_A2y=|F(P$Bg!6`uT z3lk9RwvR{&A0_Z9f*YstJ?$&D#dR4bVQ<2AWOV2%d`{(+f6pCw9%XW+vF z<|teE2SP5l5MUnVUm 0.0, alpha*w*np.exp(-w/wc), 0.0) + + def Jw2(w, wc, tau, D): + s = 0.5 + fw = None + if(D == 1): + fw = np.cos(w*tau) + if(D == 3): + fw = np.sin(w*tau)/(w*tau) + return np.power(w, D)/np.power(wc, D-1)*np.exp(-np.abs(w/wc))*(1-fw)/2.0#+np.cos(w*2*tau))/4.0 + #return w*np.exp(-np.abs(w/wc))*(1+np.cos(w*tau))/2.0#+np.cos(w*2*tau))/4.0 + #return np.pi/2.0*tau*np.power(wc, 1-s)*np.power(np.abs(w), s)*np.exp(-np.abs(w)/wc)#*(1+np.cos(w*tau))/2.0 + + def Sw(w, beta, J, *args): + return J(w, *args)*0.5*(1+1.0/np.tanh(beta*w/2.0)) + + def Sw0(w, J, *args): + return np.where(w > 0, J(w, *args), 0) + wc = 1 + g = wc*2 + + #Z1 = 1-np.logspace(-8, -2, 100) + #Z1 = np.concatenate((Z1, 1+np.logspace(-8, -2, 100))) + Z1 = np.linspace(1e-8, 20, 2000)#np.concatenate((Z1, np.linspace(1e-8, 5, 2000))) + #Z1 = + #Z1 = np.logspace(-2, np.log(1000)/np.log(10), 2000)#, np.linspace(1.01e-1, 20*wc, 1800))) + #Z1 = np.concatenate((np.logspace(-8, -1, 200), np.linspace(1.01e-1, 20*wc, 1800))) + #nZ1 = -np.concatenate((np.logspace(-8, -1, 200), np.linspace(1.01e-1, 20*wc, 100))) + nZ1 = -np.flip(Z1) + Z = np.concatenate((nZ1, Z1)) + + alpha = 1.0 + wc = 1 + + Zv = np.linspace(-20*wc, 20*wc, 100000) + func1, p, r, z = AAA_algorithm(lambda x : Jexp(x, alpha, wc), Z, nmax = 500, tol=1e-4) + print(p, r/(2.0*np.pi), z) + p=1.0j*p + Zv = np.linspace(-20*wc, 20*wc, 100000) + plt.plot(Zv, func1(Zv)-Jexp(Zv, alpha, wc)) + plt.plot(Z, func1(Z)-Jexp(Z, alpha, wc)) + #plt.plot(Zv, ) + plt.show() + + fv = func1(Zv) + f2v = Jexp(Zv, alpha, wc) + plt.plot(Zv, np.real(f2v), 'k-') + #plt.plot(Zv, np.imag(f2v), 'k--') + plt.plot(Zv, np.real(fv), 'r--') + #plt.plot(Zv, np.imag(fv), 'r--') + plt.show() + #func2 = AAA_algorithm(lambda x : Sw0(x,Jw2, 1, 10), Z, tol=1e-8,nmax=1000) + #func3 = AAA_algorithm(lambda x : Sw0(x,Jw2, 1, 10), Z, tol=1e-12,nmax=1000) + #plt.show() + + #Z = np.exp(np.linspace(-0.5,0.5+15.0j*np.pi,1000)) + #func = lambda x : np.tan(np.pi*x/2.0) + #func1, p, r, z = AAA_algorithm(func, Z, tol=1e-13,nmax=100) + #plt.show() + #print(p) + + #fig = plt.figure() + #ax = fig.add_subplot(projection='3d') + #ax.scatter(Z.real, Z.imag, np.abs(func(Z))) + #plt.show() + #exit(1) + fig = plt.figure() + ax = fig.add_subplot(projection='3d') + colors = np.where(r<0, 'r', 'b') + #ax.scatter(p.real, p.imag)#, np.log(np.abs(r)/(2.0*np.pi))/np.log(10), c=colors.ravel()) + ax.scatter(p.real, p.imag, np.log(np.abs(r))/np.log(10)) + + #pd = np.genfromtxt("poles.csv") + #poles2 = pd[:, 0] + pd[:, 1]*1.0j + #res2 = pd[:, 2] + pd[:, 3]*1.0j + #colors = np.where(res2<0, 'orange', 'cyan') + #ax.scatter(poles2.real, poles2.imag)#, np.log(np.abs(res2))/np.log(10), c=colors.ravel()) + + #inds = p.imag >= 0 + #p = p[inds] + #r = r[inds] + + for i in range(p.shape[0]): + print(p[i], r[i]/(2.0*np.pi)) + print(p.shape[0]) + + t = np.linspace(0, 100, 10000) + def Ct(p, r, t): + inds = p.real <= 0 + p = p[inds] + r = r[inds] + P, T = np.meshgrid(p, t) + return 1.0j*np.exp(P*T)@r + + plt.show() + #plt.plot(t, np.real(1.0/(t+1.0j)/(2*np.pi))) + plt.plot(t, np.real(Ct(p, r, t))) + plt.plot(t, np.imag(Ct(p, r, t))) + plt.plot(t, np.real(alpha*wc*wc/(1-1.0j*wc*t)**2)) + plt.show() + + #plt.plot(t, np.real(np.pi*alpha/2.0 * 200 *np.sqrt(np.pi)/np.power(1-20j*t, 1.5)), '--') + #plt.plot(t, np.imag(np.pi*alpha/2.0 * 200 *np.sqrt(np.pi)/np.power(1-20j*t, 1.5)), '--') + + + #plt.plot(t, np.real(Ct(poles2, res2, t))*2.0*np.pi) + #plt.plot(t, np.imag(Ct(poles2, res2, t))*2.0*np.pi) + #plt.plot(t, np.real(1.0/(8*np.pi)*wc*((-1.0+1.0j*t*wc)/(-wc*wc*tau*tau+(wc*t+1.0j)*(wc*t+1.0j) ) + 2.0j/(wc*t+1.0j)))) + + #plt.plot(t, np.real(1.0/(8*np.pi)*wc*( (-1.0+1.0j*t*wc)/(-wc*wc*tau*tau*4+(wc*t+1.0j)*(wc*t+1.0j) ) + (-1.0+1.0j*t*wc)/(-wc*wc*tau*tau+(wc*t+1.0j)*(wc*t+1.0j) ) + 2.0j/(wc*t+1.0j)))) + #plt.plot(t, np.imag(1.0/(t+1.0j)/(2*np.pi))) + #plt.plot(t, np.imag(Ct(p, r, t))) + + #plt.scatter(z.real, z.imag)# np.abs(func(Z))) + #plt.plot(Zv, func1(Zv)) + #plt.plot(Zv, func3(Zv)) + +if __name__ == "__main__": + test() diff --git a/engine/heom/bath_data.py b/engine/heom/bath_data.py new file mode 100644 index 0000000..6869dfc --- /dev/null +++ b/engine/heom/bath_data.py @@ -0,0 +1,102 @@ +import numpy as np +from .aaa import AAA_algorithm + +def softmspace(start, stop, N, beta = 1, endpoint = True): + start = np.log(np.exp(beta*start)-1)/beta + stop = np.log(np.exp(beta*stop)-1)/beta + + dx = (stop-start)/N + if(endpoint): + dx = (stop-start)/(N-1) + + return np.log(np.exp(beta*(np.arange(N)*dx + start))+1)/beta + +def generate_grid_points(N, wc, wmin=1e-9): + Z1 = softmspace(wmin, 20*wc, N) + nZ1 = -np.flip(Z1) + Z = np.concatenate((nZ1, Z1)) + return Z + + +def AAA_to_HEOM(p, r): + pp = p*1.0j + rr = -1.0j*r/(np.pi) + inds = pp.real > 0 + pp = pp[inds] + rr = rr[inds] + return rr, pp + +def setup_heom_correlation_functions(Sw, Z1, nmax = 500, aaa_tol = 1e-4): + #first compute the aaa decomposition of the spectral function + func1, p, r, z = AAA_algorithm(Sw, Z1, nmax=nmax, tol=aaa_tol) + + #and convert that to the heom correlation function coefficients + dk, zk = AAA_to_HEOM(p, r) + + #return the function for optional plotting as well as the coefficients + return func1, dk, zk + +class heom_bath: + def __init__(self, Scoup, Sw, L, Lmin = None, aaa_support_points = None, wmax = None, wmin = None, aaa_tol = 1e-3, Naaa=1000, aaa_nmax=500, scale_factor=None): + self.Scoup = Scoup + self.Sw = Sw + self.L = L + self.Lmin = None + self._aaa_support_points = aaa_support_points + self.wmax = wmax + self.wmin = wmin + self.aaa_tol = aaa_tol + self.Naaa = Naaa + self.aaa_nmax=500 + self.scale_factor = scale_factor + + def aaa_support_points(self): + if(self._aaa_support_points == None): + wmax = self.wmax + if(self.wmax == None): + wmax = 1 + + wmin = self.wmin + if(self.wmin == None): + wmin = 1e-8 + + return generate_grid_points(self.Naaa, wmax, wmin=wmin) + + elif isinstance(self._aaa_support_points, (list, np.ndarray)): + if isinstance(self._aaa_support_points, list): + return np.array(self._aaa_support_points) + else: + return self._aaa_support_points + + + def discretise(self, output_fitting=False): + Z1 = self.aaa_support_points() + Sw_aaa, dk, zk = setup_heom_correlation_functions(self.Sw, Z1, nmax=self.aaa_nmax, aaa_tol=self.aaa_tol) + + bath_dict = { + "S" : self.Scoup, + "d" : dk, + "z" : zk, + "L" : self.L + } + + if self.Lmin is not None: + bath_dict["Lmin"]=self.Lmin + + if self.scale_factor is not None: + bath_dict["sf"]=self.scale_factor + + if output_fitting: + wmax = self.wmax + if(self.wmax == None): + wmax = 1 + Zv = np.linspace(-20*wmax, 20*wmax, 100000) + Swa = Sw_aaa(Zv) + Swb = self.Sw(Zv) + + bath_dict["wf"] = Zv + bath_dict["Sw"] = Swb + bath_dict["Sw_fit"] = Swa + + return bath_dict + diff --git a/engine/heom/heomif.py b/engine/heom/heomif.py new file mode 100644 index 0000000..1c16373 --- /dev/null +++ b/engine/heom/heomif.py @@ -0,0 +1,175 @@ +from .aaa import AAA_algorithm +from .mps import mps + +import numpy as np +import scipy as sp + + +def commutator(L): + return np.kron(L, np.identity(L.shape[0])) - np.kron(np.identity(L.shape[0]), L.T) + +def anti_commutator(L): + return np.kron(L, np.identity(L.shape[0])) + np.kron(np.identity(L.shape[0]), L.T) + + +def Lkp(nbose, dk, S, s=0.5): + b1 = np.zeros((nbose, nbose), dtype=np.complex128) + + for i in range(nbose-1): + b1[i, i+1] = np.sqrt((i+1.0))*np.sqrt(np.abs(dk)) + + Scomm = commutator(S) + return np.kron(Scomm, b1) + + +def Lkm(nbose, dk, S, mind = True, s=0.5): + b1 = np.zeros((nbose, nbose), dtype=np.complex128) + + coeff = 1 + if not mind: + coeff = -1.0 + for i in range(nbose-1): + b1[i+1, i] = coeff*np.sqrt((i+1.0))*dk/np.sqrt(np.abs(dk)) + + Sop = None + if mind: + Sop = np.kron(S, np.identity(S.shape[0])) + else: + Sop = np.kron(np.identity(S.shape[0]), S.T) + return np.kron(Sop, b1) + +def nop(D2, nbose, zk): + return -1.0j*np.kron(np.identity(D2), np.diag((np.arange(nbose))*zk)) + + +def mode_operator(nbose, dk, zk, S, mind): + op = None + s = 0.5 + if mind: + op = Lkp(nbose, dk, S, s=s) + Lkm(nbose, dk, S, mind=mind, s= s) + nop(S.shape[0]**2, nbose, zk) + else: + op = Lkp(nbose, np.conj(dk), S, s=s) + Lkm(nbose, np.conj(dk), S, mind=mind, s=s) + nop(S.shape[0]**2, nbose, np.conj(zk)) + return op + + +def Mk(nbose, dk, zk, S, mind, dt, nf=None): + op = None + if(nf > nbose): + op = mode_operator(nf, dk, zk, S, mind) + s = S.shape[0]*S.shape[1] + expm = sp.linalg.expm(-1.0j*dt*op).reshape(s, nf, s, nf) + return expm[:, :nbose, :, :nbose].reshape(s*nbose, s*nbose) + else: + op = mode_operator(nbose, dk, zk, S, mind) + return sp.linalg.expm(-1.0j*dt*op) + +def compute_dimensions(S, dk, zk, L, Lmin = None): + ds = np.ones(2*len(dk)+1, dtype = int) + ds[0] = S.shape[0]*S.shape[1] + + minzk = np.amin(np.real(zk)) + if(Lmin == None): + for i in range(len(dk)): + nb = L + ds[2*i+1] = nb + ds[2*i+2] = nb + else: + for i in range(len(dk)): + nb = max(int(L*minzk/np.real(zk[i])), Lmin) + ds[2*i+1] = nb + ds[2*i+2] = nb + + return ds + +#build the short time HEOM propagator +def build_propagator_matrices(S, dk, zk, dt, L, Lmin=None, sf = 1): + ds = compute_dimensions(S, dk, zk, L, Lmin=Lmin) + + Uks = [] + for i in range(len(dk)): + nb = ds[2*i+1] + Uks.append(Mk(nb, dk[i], zk[i], S, True, dt/2.0, nf = sf*nb)) + Uks.append(Mk(nb, dk[i], zk[i], S, False, dt/2.0, nf = sf*nb)) + return Uks + + +def HEOM_bath_propagator(bath, dt): + Lmin = None + sf = 1 + if 'Lmin' in bath.keys(): + Lmin = bath['Lmin'] + if 'sf' in bath.keys(): + sf = bath['sf'] + return build_propagator_matrices(bath['S'], bath['d'], bath['z'], dt, bath['L'], Lmin=Lmin, sf=sf) + +def HEOM_propagator(baths, dt): + #if we have a list of baths + if isinstance(baths, list): + Uks = [] + for bath in baths: + Uks = Uks + HEOM_bath_propagator(bath, dt) + return Uks + + elif isinstance(baths, dict): + return HEOM_bath_propagator(baths, dt) + + +def apply_propagator(Us, Uks, A, method='naive', tol=None, nbond = None): + #apply non local two site gates. Here we perform swap operations as we go, unless we are applying the last gate + for i, Uk in enumerate(Uks): + if(i+1 == len(Uks)): + A.apply_two_site(Uk, i, i+1, method=method, tol=tol, nbond=nbond) + else: + A.apply_bond_tensor_and_swap(Uk, i, dir = 'right', tol=tol, nbond=nbond) + + A.apply_one_site(Us, -2) + + for i, Uk in reversed(list(enumerate(Uks))): + if(i+1 == len(Uks)): + A.apply_two_site(Uk, i, i+1, method=method, tol=tol, nbond=nbond) + else: + A.apply_bond_tensor_and_swap(Uk, i, dir = 'left', tol=tol, nbond=nbond) + + return A + + + +def setup_HEOM_ados(rho0, prop): + nliouv = rho0.flatten().shape[0] + Nmodes = len(prop)+1 + d = np.zeros(Nmodes, dtype = int) + for i, Uk in enumerate(prop): + d[i+1] = Uk.shape[0]//nliouv + d[0] = nliouv + + #setup the ado mps + A = mps(chi = np.ones(len(prop), dtype=int), d=d, init='zeros', dtype=np.complex128) + + A[0][0, :, 0] = rho0.flatten() + for i in range(1, len(A)): + A[i][0, 0, 0] = 1.0 + A.orthogonalise() + return A + + +def extract_rho(A): + T = None + for i in reversed(range(1, len(A))): + if not isinstance(T, np.ndarray): + T = A[i][:, 0, 0] + else: + M = A[i][:, 0, :] + T = M@T + + Mi = A[0][:, :, :] + T = np.tensordot(Mi, T, axes=([2], [0])) + return T[0, :] + + +def Cr(func, t, tol=1e-8, limit=1000, workers=1, dx = 1e-12): + return sp.integrate.quad_vec(lambda w : np.real(func(w)*np.exp(-1.0j*t*w)), dx, np.inf, limit=1000)[0]/np.pi + sp.integrate.quad_vec(lambda w : np.real(func(w)*np.exp(-1.0j*t*w)), -np.inf, -dx, limit=1000)[0]/np.pi + sp.integrate.quad_vec(lambda w : np.real(func(w)*np.exp(-1.0j*t*w)), -dx, dx, limit=1000, points=[0])[0]/np.pi + +def Ci(func, t, tol=1e-8, limit=1000, workers=1, dx = 1e-12): + return sp.integrate.quad_vec(lambda w : np.imag(func(w)*np.exp(-1.0j*t*w)), dx, np.inf, limit=1000)[0]/np.pi + sp.integrate.quad_vec(lambda w : np.imag(func(w)*np.exp(-1.0j*t*w)), -np.inf, -dx, limit=1000)[0]/np.pi + sp.integrate.quad_vec(lambda w : np.imag(func(w)*np.exp(-1.0j*t*w)), -dx, dx, limit=1000, points=[0])[0]/np.pi + + diff --git a/engine/heom/mps.py b/engine/heom/mps.py new file mode 100644 index 0000000..82de648 --- /dev/null +++ b/engine/heom/mps.py @@ -0,0 +1,780 @@ +import numpy as np +from .mps_utils import * + + +#a basic shared memory mps implementation. This contains a number +class mps: + #construct from array of bond dimensions and local hilbert space dimensions + def __init__(self, chi = None, d = None, N = None, chil : int = 1, chir : int = 1, dtype = np.double, init = 'zeros', do_orthog=True): + self.dtype = dtype + if isinstance(chi, (np.ndarray, list)) and isinstance(d, (np.ndarray, list)): + self.build_from_arrays(chi, d, chil = chil, chir = chir, dtype = dtype, init = init, do_orthog=do_orthog) + elif isinstance(chi, int) and isinstance(d, (np.ndarray, list)): + N = len(d) + chi_arr = np.ones(N-1, dtype=int)*chi + self.build_from_arrays(chi_arr, d, chil = chil, chir = chir, dtype = dtype, init = init, do_orthog=do_orthog) + elif isinstance(chi, (np.ndarray, list)) and isinstance(d, int): + N = len(chi)+1 + d_arr = np.ones(N, dtype=int)*d + self.build_from_arrays(chi, d_arr, chil = chil, chir = chir, dtype = dtype, init = init, do_orthog=do_orthog) + else: + if N == None: + self._tensors = None + self._is_valid = False + self._boundary_vects = None + self._orth_centre = None + else: + chi_arr = np.ones(N-1, dtype=int)*chi + d_arr = np.ones(N, dtype=int)*d + self.build_from_arrays(chi_arr, d_arr, chil = chil, chir = chir, dtype = dtype, init = init, do_orthog=do_orthog) + + self._is_conjugated = False + + def build_from_arrays(self, chi, d, chil = 1, chir = 1, dtype = np.double, init = 'zeros', do_orthog=True): + Nchi = None + Nd = None + if not isinstance(chi, (list, np.ndarray)): + raise TypeError("Bond dimension variable is an invalid type.") + if not isinstance(d, (list, np.ndarray)): + raise TypeError("Local Hilbert Space Dimension variable is an invalid type.") + + Nchi = len(chi) + Nd = len(d) + + if not (Nchi + 1 == Nd): + raise RuntimeError("bond dimension and local hilbert space arrays are not compatible.") + + self._tensors = [None]*Nd + if Nd > 1: + self._tensors[0] = np.zeros((chil, d[0], chi[0]), dtype=dtype) + for i in range(1, Nd-1): + self._tensors[i] = np.zeros((chi[i-1], d[i], chi[i]), dtype=dtype) + self._tensors[-1] = np.zeros((chi[-1], d[-1], chir), dtype=dtype) + elif Nd == 1: + self._tensors[0] = np.zeros((chil, d[0], chir), dtype=dtype) + + self._orth_centre = None + self._is_valid = True + self._boundary_vects = None + self.init_values(dtype=dtype, init=init, do_orthog=do_orthog) + + def init_values(self, dtype = np.double, init='zeros', do_orthog=True): + if isinstance(init, str): + if init == 'zeros': + for A in self._tensors: + A = np.zeros(A.shape, dtype=dtype) + + elif init == 'random': + for A in reversed(self._tensors): + if(dtype == np.complex128): + A += np.random.normal(0, 1, size = (np.prod(A.shape), 2)).view(np.complex128).flatten().reshape(A.shape) + else: + A += np.random.normal(0, 1, size=A.shape) + + if(do_orthog): + for i in reversed(range(self.nbonds())): + self.shift_bond(i, dir='left') + norm = np.dot(np.conj(self._tensors[i].flatten()), self._tensors[i].flatten()) + self._tensors[i] /= np.sqrt(norm) + + self._orth_centre = 0 + + elif isinstance(init, (np.ndarray, list)): + if(len(init) != self.nsites()): + raise ValueError("Invalid init array") + + for i in range(self.nsites()): + self._tensors[i] = np.zeros(self._tensors[i].shape, dtype=dtype) + ind = int(init[i]) + if(ind >= self._tensors[i].shape[1] or ind < 0): + raise IndexError("Index out of bounds for state vector initialisation") + + self._tensors[i][0, ind, 0] = 1.0 + if(do_orthog): + self.orthogonalise() + + + #access tensor by index + def __getitem__(self, i): + if isinstance(i, slice): + return [self[ii] for ii in range(*i.indices(len(self)))] + elif isinstance(i, int): + i = mapint(i, len(self)) + if self._is_conjugated: + return np.conj(self._tensors[i]) + else: + return self._tensors[i] + else: + raise TypeError("Invalid argument type") + return None + + def __getslice__(self, i): + return self.__getitem__(i) + + #set tensor by index + def __setitem__(self, i, v): + #first we check that v is a valid type + if not isinstance(v, np.ndarray): + raise TypeError("Invalid type for setting item.") + if not v.ndim == 3: + raise TypeError("Invalid type for setting item.") + + if isinstance(i, slice): + for ii in range(*i.indices(len(self))): + self._tensors[ii] = v + elif isinstance(i, int): + i = mapint(i, len(self)) + self._tensors[i] = v + else: + raise TypeError("Invalid argument type") + + self._orth_centre = None + self._is_valid = None + self.is_valid() + + def __setslice__(self, i, v): + self.__setitem__(i, v) + + def is_ortho(self): + return not self._orth_centre == None + + def __len__(self): + return len(self._tensors) + + def nsites(self): + return len(self) + + def nbonds(self): + return len(self)-1 + + def sweep_order(self): + return list(range(len(self))) + + def conj(self): + ret = mps() + ret._tensors = self._tensors + ret._orth_centre = self._orth_centre + ret._boundary_vects = self._boundary_vects + ret._is_valid = self._is_valid + ret._is_conjugated = not self._is_conjugated + return ret + + def isconjugated(self): + return self._is_conjugated + + def shape(self, i): + if isinstance(i, slice): + return [self.shape(ii) for ii in range(*i.indices(len(self)))] + elif isinstance(i, int): + i = mapint(i, len(self)) + return self._tensors[i].shape + else: + raise TypeError("Invalid argument type") + + def local_dimension(self, i): + if isinstance(i, slice): + return [self.local_dimension(ii) for ii in range(*i.indices(len(self)))] + elif isinstance(i, int): + i = mapint(i, len(self)) + return self._tensors[i].shape[1] + else: + raise TypeError("Invalid argument type") + + def expand_local_dimension(self, i, d): + if isinstance(i, int): + i = mapint(i, self.nsites()) + if(d > self._tensors[i].shape[1]): + npad = d - self._tensors[i].shape[1] + self._tensors[i] = np.pad(self._tensors[i], ((0,0), (0, npad), (0,0) ), 'constant', constant_values=(0)) + else: + raise TypeError("Invalid argument type") + + def bond_dimension(self, i): + if isinstance(i, slice): + return [self.bond_dimension(ii) for ii in range(*i.indices(len(self)))] + elif isinstance(i, int): + i = mapint(i, self.nbonds()) + return self._tensors[i].shape[2] + else: + raise TypeError("Invalid argument type") + + def maximum_bond_dimension(self): + bd = 0 + for i in range(self.nbonds()): + if self._tensors[i].shape[2] > bd: + bd = self._tensors[i].shape[2] + return bd + + def expand_bond(self, i, chi): + if isinstance(i, int): + i = mapint(i, self.nbonds()) + if not self.shape(i)[2] == self.shape(i+1)[0]: + raise RuntimeError("Cannot expand bond, the specified bond is invalid.") + if(chi > self.shape(i)[2]): + npad = chi - self.shape(i)[2] + self._tensors[i] = np.pad(self._tensors[i], ((0,0), (0, 0), (0,npad) ), 'constant', constant_values=(0)) + self._tensors[i+1] = np.pad(self._tensors[i+1], ((0,npad), (0, 0), (0,0) ), 'constant', constant_values=(0)) + else: + raise TypeError("Invalid argument type") + + #a function for taking an invalid object and sanitising the bonds to make sure that the object is a valid size. + def sanitise(self): + #if we are valid then we don't need to do anything. + if not self.is_valid(): + for bi in range(self.nbonds()): + chil = self.shape(bi)[2] + chir = self.shape(bi+1)[0] + if not chil == chir: + chi = max(chil, chir) + if(chi > chil): + npad = chi - chil + self._tensors[bi] = np.pad(self._tensors[i], ((0,0), (0, 0), (0,npad) ), 'constant', constant_values=(0)) + else: + npad = chi - chir + self._tensors[i+1] = np.pad(self._tensors[i+1], ((0,npad), (0, 0), (0,0) ), 'constant', constant_values=(0)) + self._is_valid = True + + + #function for ensuring that the mps is currently valid + def is_valid(self): + if not self._is_valid == None: + return self._is_valid + else: + iv = True + for i in range(len(self._tensors)): + if not self.shape(i-1)[2] == self.shape(i)[0]: + iv = False + self._is_valid = iv + return self._is_valid + + def orthogonalise(self): + for i in reversed(range(self.nbonds())): + self.shift_bond(i, dir='left') + self._orth_centre = 0 + + def truncate(self, tol = None, nbond=None): + self.shift_orthogonality(-1) + self.shift_orthogonality(0, tol=tol, nbond=nbond) + + def shift_orthogonality(self, i, tol = None, nbond = None): + i = mapint(i, self.nsites()) + if self._orth_centre == None: + self.orthogonalise() + + if i < self._orth_centre: + for bi in reversed(range(i, self._orth_centre)): + self.shift_left(tol = tol, nbond = nbond) + + elif i > self._orth_centre: + for bi in range(self._orth_centre, i): + self.shift_right(tol = tol, nbond = nbond) + + if(i != self._orth_centre): + raise RuntimError("the orthogonality centre has not been shifted to the correct position.") + + + def normalise(self): + norm = None + if(self._orth_centre != None): + oc = self._orth_centre + norm = np.dot(np.conj(self._tensors[oc].flatten()), self._tensors[oc].flatten()) + self._tensors[oc] /= np.sqrt(norm) + else: + norm = contract(self, self) + self._tensors[0] /= np.sqrt(norm) + return norm + + + def expand(self, T, boundary, is_orthogonalised = False): + if not isinstance(T, np.ndarray): + raise RuntimeError("Cannot expand mps as input tensor is invalid.") + if T.ndim != 3: + raise RuntimeError("Input tensors is not the correct dimension.") + + if(boundary == 'left'): + if(T.shape[2] != self._tensors[0].shape[0]): + raise RuntimeError("Cannot expand MPS boundary tensor has incompatible shape.") + self._tensors.insert(0, T) + + if not is_orthogonalised: + self._orth_centre = None + else: + self._orth_centre += 1 + + elif(boundary == 'right'): + if(T.shape[0] != self._tensors[self.nbonds()].shape[2]): + raise RuntimeError("Cannot expand MPS boundary tensor has incompatible shape.") + + self._tensors.append(T) + else: + raise RuntimeError("Failed to expand MPS. Boundary type not recognised") + + def pop(self, boundary): + if(boundary == 'left'): + if not self._orth_centre == None: + if(self._orth_centre > 0): + self._orth_centre -= 1 + return self._tensors.pop(0) + + elif(boundary == 'right'): + if not self._orth_centre == None: + if(self._orth_centre == self.nbonds() ): + self._orth_centre -= 1 + return self._tensors.pop(-1) + + else: + raise RuntimeError("Failed to pop boundary tensor from MPS. Boundary type not recognised") + + + def shift_left(self, tol = None, nbond = None): + if not self.is_valid(): + raise RuntimeError("The object does not represent a valid MPS. Unable to perform transformation operations to the MPS.") + if self._orth_centre == None: + raise RuntimeError("The object does not have an orthogonality centre to shift") + if self._orth_centre == 0: + raise RuntimeError("Orthogonality Centre cannot be shifted left") + self.shift_bond(self._orth_centre - 1, dir='left', tol = tol, nbond = nbond) + + def shift_right(self, tol = None, nbond = None): + if not self.is_valid(): + raise RuntimeError("The object does not represent a valid MPS. Unable to perform transformation operations to the MPS.") + if self._orth_centre == None: + raise RuntimeError("The object does not have an orthogonality centre to shift") + if self._orth_centre == self.nbonds(): + raise RuntimeError("Orthogonality Centre cannot be shifted left") + self.shift_bond(self._orth_centre , dir='right', tol = tol, nbond = nbond) + + + def shift(self, dir, tol = None, nbond = None): + if dir == 'left': + self.shift_left(tol=tol, nbond=nbond) + elif dir == 'right': + self.shift_right(tol=tol, nbond=nbond) + else: + raise RuntimeError("Failed to shift bond incorrect direction.") + + #updates the MPS so that the site tensors are isometries and return the non-orthogonal bond matrix + def schmidt_decomposition(self, dir, tol=None, nbond=None, chimin=None): + if self._orth_centre == None: + raise RuntimeError("The schmidt decomposition function requires the MPS to be in a mixed canonical form.") + if (self._orth_centre == 0 and dir == 'left') or (self._orth_centre + 1 == len(self) and dir == 'right'): + raise RuntimeError("Unable to perform specified decomposition we are at the bounds of the MPS.") + + bind = None + if dir == 'left': + bind = self._orth_centre-1 + else: + bind = self._orth_centre + + il = bind + ir = bind+1 + + self._tensors[il], R, self._tensors[ir] = local_canonical_form(self._tensors[il], self._tensors[ir], dir, il, self._orth_centre, tol = tol, nbond = nbond, chimin = chimin) + return R + + def shift_bond(self, bind, dir='right', tol=None, nbond=None, chimin=None): + if not self.is_valid(): + raise RuntimeError("The object does not represent a valid MPS. Unable to perform transformation operations to the MPS.") + + bind = mapint(bind, self.nbonds()) + + #get the indices of the two sites that will be modified by the operation + il = bind + ir = bind+1 + + self._tensors[il], self._tensors[ir] = shift_mps_bond(self._tensors[il], self._tensors[ir], dir, il, self._orth_centre, tol = tol, nbond = nbond, chimin = chimin) + self._orth_centre = update_mps_ortho_centre(il, ir, self._orth_centre, dir) + + def state_internal(self, init): + if isinstance(init, (np.ndarray, list)): + if(len(init) != self.nsites()): + raise ValueError("Invalid init array") + M = None + for i in range(self.nsites()): + ind = int(init[i]) + if(ind >= self._tensors[i].shape[1] or ind < 0): + raise IndexError("Index out of bounds for state vector initialisation") + if not isinstance(M, np.ndarray): + M = self._tensors[i][:, ind, :] + else: + M = M @ self._tensors[i][:, ind, :] + return M + else: + raise RuntimeError("Not enough components for state") + + def state(self, init): + M = self.state_internal(init) + if(M.shape[0] == 1 and M.shape[1] == 1): + return M[0, 0] + else: + return M + + def __str__(self): + if(self._orth_centre != None): + return 'MPS: tensors: %s \n orth centre: %d'%(self._tensors, self._orth_centre) + else: + return 'MPS: tensors: %s'%(self._tensors) + + def __imul__(self, x): + if(self._orth_centre != None): + self._tensors[self._orth_centre] *= x + else: + self._tensors[0] *= x + return self + + def __itruediv__(self, x): + if(self._orth_centre != None): + self._tensors[self._orth_centre] /= x + else: + self._tensors[0] /= x + return self + + #function for applying a general one-site operator to the MPS + def apply_one_site(self, M, i, shift_orthogonality = False): + if self._is_conjugated: + M = np.conj(M) + i = mapint(i, self.nsites()) + + dM = M.shape[1] + dims = self._tensors[i].shape + if(dM != dims[1]): + raise RuntimeError("The one site operator and MPS site tensor do not have compatible dimensions.") + + if(self._orth_centre != None and shift_orthogonality): + self.shift_orthogonality(i) + + self._tensors[i] = np.transpose(np.tensordot(M, self._tensors[i], axes=([1], [1])), axes=[1, 0, 2]) + + if(self._orth_centre != i): + self._orth_centre = None + #self._tensors[i] = np.einsum('ij, ajb -> aib', M, self._tensors[i]) + + #function for applying a two site tensor across a bond and swapping the resultant output indices of the bond + def apply_bond_tensor_and_swap(self, M, bi, dir='right', tol=None, nbond=None, optol=None): + bi = mapint(bi, self.nbonds()) + + mi = None + #if we are going the 'right' direction then we have bi as the first index and bi+1 as the second index + if dir == 'right': + mi = [bi, bi+1] + else: + mi = [bi+1, bi] + + #now we have set up the state in a suitable form for constructing the + if(self._orth_centre == None): + self.orthogonalise() + + self.shift_orthogonality(bi) + dims = [self.local_dimension(m) for m in mi] + + #now we convert the n-site object into an MPO + Mt, inds, ds = permute_nsite_dims(M, mi, dims) + Mp = nsite_mpo(Mt, ds, order='clockwise', tol=optol) + + i = bi + j = bi+1 + + c=0 + for ind in range(i, j+1): + self._tensors[ind] = mps.apply_MPO_node(Mp[c], self._tensors[ind], order='clockwise') + c = c+1 + + Mv = np.transpose(self.contract_bond(bi), axes=[0, 2, 1, 3]) + self.decompose_two_site(Mv, bi, dir=dir, tol=tol, nbond=nbond) + if dir == 'right': + self._orth_centre = bi+1 + else: + self._orth_centre = bi + + #performs a swap operation of physical indices across a bond + def swap_bond_indices(self, bi, dir='right', tol=None, nbond=None): + #ensure that the orthogonality centre is at site bi (the left most of the two sites involved in the bond) + bi = mapint(bi, self.nbonds()) + self.shift_orthogonality(bi) + M = np.transpose(self.contract_bond(bi), axes=[0, 2, 1, 3]) + self.decompose_two_site(M, bi, dir=dir, tol=tol, nbond=nbond) + if dir == 'right': + self._orth_centre = bi+1 + else: + self._orth_centre = bi + + #contracts an MPS across a bond to the two-site tensor object + def contract_bond(self, bi): + bi = mapint(bi, self.nbonds()) + return np.tensordot(self._tensors[bi], self._tensors[bi+1], axes=([2], [0])) + + def decompose_two_site(self, M, bi, dir = 'right', tol = None, nbond = None): + bi = mapint(bi, self.nbonds()) + dims = M.shape + + il = bi + ir = bi+1 + + Mm = M.reshape((dims[0]*dims[1], dims[2]*dims[3])) + + Q, S, Vh = np.linalg.svd(Mm, full_matrices=False, compute_uv=True) + nsvd = determine_truncation(S, tol = tol, nbond = nbond) + + Q = Q[:, :nsvd] + S = np.diag(S[:nsvd]) + Vh = Vh[:nsvd, :] + + if dir == 'right': + R = S @ Vh + + self._tensors[il] = Q.reshape( (dims[0], dims[1], nsvd)) + self._tensors[ir] = R.reshape( (nsvd, dims[2], dims[3])) + + if self._orth_centre == il: + self._orth_centre = ir + + elif dir == 'left': + R = Q @ S + self._tensors[il] = R.reshape( (dims[0], dims[1], nsvd)) + self._tensors[ir] = Vh.reshape( (nsvd, dims[2], dims[3])) + + if self._orth_centre == ir: + self._orth_centre = il + else: + ValueError("Invalid dir argument") + + def apply_MPO_node(Mt, nt, order = 'mpo'): + ret = None + if(order == 'mpo'): + ret = np.transpose(np.tensordot(Mt, nt, axes=([2], [1])), (0, 3, 1, 2, 4)) + #ret = np.einsum('aijb, mjn -> amibn', Mt, nt) + else: + ret = np.transpose(np.tensordot(Mt, nt, axes=([3], [1])), (0, 3, 1, 2, 4)) + #ret = np.einsum('aibj, mjn -> amibn', Mt, nt) + + d = ret.shape + return (ret.reshape((d[0]*d[1], d[2], d[3]*d[4]))) + + #function for applying a general two-site operator to the MPS + def apply_two_site(self, M, i, j, method='naive', tol = None, nbond = None, optol=None): + self._apply_nsite(M, [i, j], method=method, tol=tol, nbond=nbond, optol=optol) + + def _apply_nsite(self, M, inds, method='naive', tol = None, nbond = None, optol = None): + mi = [mapint(x, len(self)) for x in inds] + + if len(mi) > len(self): + raise ValueError("Index array too large.") + if len(mi) != len(set(mi)): + raise ValueError("Index array contained duplicates.") + + if(M.shape[0] != M.shape[1]): + raise ValueError("Input operator is incorrect size.") + + mdim = 1 + dims = [self.local_dimension(m) for m in mi] + for i, m in enumerate(mi): + mdim = mdim * self.local_dimension(m) + + if(mdim != M.shape[0]): + raise ValueError("Input operator and inds array are incompatible for this MPS.") + + if self._is_conjugated: + M = np.conj(M) + + #now we have set up the state in a suitable form for constructing the + if(self._orth_centre == None): + self.orthogonalise() + + #now we convert the n-site object into an MPO + Mt, inds, ds = permute_nsite_dims(M, mi, dims) + Mp = nsite_mpo(Mt, ds, order='clockwise', tol=optol) + + i = inds[0] + j = inds[-1] + + c = 0 + + if method == "naive": + nbd = 1 + for ind in range(i, j+1): + if(ind == inds[c]): + self._tensors[ind] = mps.apply_MPO_node(Mp[c], self._tensors[ind], order='clockwise') + nbd = Mp[c].shape[2] + c = c+1 + else: + d = self._tensors[ind].shape + self._tensors[ind] = mps.apply_MPO_node(identity_pad_mpo(nbd, d[1], order='clockwise'), self._tensors[ind], order='clockwise') + + self.shift_orthogonality(j) + + #and shift back to the original orthogonality centre truncating all bonds + self.shift_orthogonality(i, tol=tol, nbond=nbond) + + else: + raise ValueError("Invalid two site mpo mps contraction scheme.") + + + #function for applying a general two-site operator to the MPS + def apply_nsite(self, M, inds, method='naive', tol = None, nbond = None, optol=None): + if isinstance(inds, int): + self.apply_one_site(M, inds) + elif isinstance(inds, list): + if len(inds) == 1: + self.apply_one_site(M, inds[0]) + else: + self._apply_nsite(M, inds, method=method, tol=tol, nbond=nbond, optol=optol) + else: + raise RuntimeError("Invalid index object passed to apply_nsite.") + + def apply_operator(self, op, method="naive", tol=None, nbond=None, optol=None): + if isinstance(op, dict): + if not "op" in op: + raise ValueError("Invalid dictionary for applying operator") + else: + if "mode" in op: + self.apply_one_site(op["op"], op["mode"]) + elif "modes" in op: + if len(op["modes"]) == 1: + self.apply_one_site(op["op"], op["modes"][0]) + elif len(op["modes"]) > 1: + if "tol" in op: + tol = op["tol"] + if "nbond" in op: + nbond = op["nbond"] + if "method" in op: + method = op["method"] + self.apply_nsite(op["op"], op["modes"], method=method, tol=tol, nbond=nbond, optol=optol) + else: + raise ValueError("Two body or higher operators only implemented through MPO") + else: + raise ValueError("Failed to read information about modes") + + def __add__(self, other): + check_compatible(self, other, mps, mps) + + if(self.shape(0)[0] != other.shape(0)[0] or self.shape(-1)[2] != other.shape(-1)[2]): + raise ValueError("Unable to add mps objects with different exterior bond dimensions.") + chil = self.shape(0)[0] + chir = self.shape(-1)[2] + + bslice = slice(0, len(self)-1) + chis = [x+y for x,y in zip(self.bond_dimension(bslice), other.bond_dimension(bslice))] + ds = self.local_dimension(slice(0, len(self))) + + ret = mps(chi=chis, d=ds, chil=chil, chir=chir, dtype = self.dtype, do_orthog=False) + + cs = self.bond_dimension(0) + ret[0][:, :, 0:cs] = self[0] + ret[0][:, :, cs:chis[0]] = other[0] + #for all interior tensors we set up the matrices as required + for i in range(1, len(self)-1): + cs = self.bond_dimension(i) + ret[i][0:cs, :, 0:cs] = self[i] + ret[i][cs:chis[i], :, cs:chis[i]] = other[i] + + cs = self.bond_dimension(-1) + ret[-1][0:cs, :, :] = self[-1] + ret[-1][cs:chis[-1], :, :] = other[-1] + + return ret + + def __sub__(self, other): + check_compatible(self, other, mps, mps) + + if(self.shape(0)[0] != other.shape(0)[0] or self.shape(-1)[2] != other.shape(-1)[2]): + raise ValueError("Unable to add mps objects with different exterior bond dimensions.") + chil = self.shape(0)[0] + chir = self.shape(-1)[2] + + bslice = slice(0, len(self)-1) + chis = [x+y for x,y in zip(self.bond_dimension(bslice), other.bond_dimension(bslice))] + ds = self.local_dimension(slice(0, len(self))) + + ret = mps(chi=chis, d=ds, chil=chil, chir=chir, dtype = self.dtype, do_orthog=False) + + cs = self.bond_dimension(0) + ret[0][:, :, 0:cs] = self[0] + ret[0][:, :, cs:chis[0]] = -other[0] + + #for all interior tensors we set up the matrices as required + for i in range(1, len(self)-1): + cs = self.bond_dimension(i) + ret[i][0:cs, :, 0:cs] = self[i] + ret[i][cs:chis[i], :, cs:chis[i]] = other[i] + + cs = self.bond_dimension(-1) + ret[-1][0:cs, :, :] = self[-1] + ret[-1][cs:chis[-1], :, :] = other[-1] + + return ret + + def overlap(A, B): + check_compatible(A, B, mps, mps); + + #now that we have checked that the sizes are correct we can perform the contraction of the MPS + res = None + for i in range(A.nsites()): + Ai = A[i] + Bi = B[i] + #contract the first node to form the res object + if(i == 0): + #res = np.einsum('ajb, cjd -> acbd', A[i], B[i]) + res = np.transpose(np.tensordot(Ai, Bi, axes=([1], [1])), axes=[0, 2, 1, 3]) + else: + #temp = np.einsum('acbd, bjm -> acmdj', res, A[i]) + temp = np.transpose(np.tensordot(res, Ai, axes=([2], [0])), axes=[0, 1, 4, 2, 3]) + res = np.tensordot(temp, Bi, axes=([3, 4], [0, 1])) + #res = np.einsum('acmdj, djn -> acmn', temp, B[i]) + if(np.prod(res.shape)== 1): + return res[0, 0, 0, 0] + else: + return res + + def matrix_element(A, M, B): + check_compatible(A, B, mps, mps); + check_compatible(A, M, mps, mpo) + + #now that we have checked that the sizes are correct we can perform the contraction of the MPS + res = None + for i in range(A.nsites()): + #contract the first node to form the res object + if(i == 0): + temp = np.transpose(np.tensordot(A[i], M[i], axes=([1], [1])), axes=(0, 2, 3, 1, 4)) + #temp = np.einsum('ajb, cjkd -> ackbd', A[i], M[i]) + res = np.transpose(np.tensordot(temp, B[i], axes=([2], [1])), axes=(0, 1, 4, 2, 3, 5)) + #res = np.einsum('ackbd, ekf -> acebdf', temp, B[i]) + else: + temp = np.einsum('acebdf, bjm -> acemjdf', res, A[i]) + temp2 = np.einsum('acemjdf, djkn -> acemnfk', temp, M[i]) + res = np.einsum('acemnfk, fko -> acemno', temp2, B[i]) + + if(np.prod(res.shape)== 1): + return res[0, 0, 0, 0, 0, 0] + else: + return res + + def contract(A, B=None, MPO=None): + if B is None: + if MPO is None: + if A.is_ortho(): + Ao = A[A._orth_centre] + return np.tensordot(np.conj(Ao), Ao, axes=([0,1,2], [0, 1, 2])) + #return np.einsum('ijk, ijk', np.conj(Ao), Ao) + else: + return mps.overlap(A.conj(), A) + else: + if isinstance(MPO, mpo): + return mps.matrix_element(A.conj(), MPO, A) + #allow for the evaluation of one site matrix elements through the contract function. + #elif isinstance(MPO, dict): + else: + if MPO is None: + return mps.overlap(B, A) + else: + return mps.matrix_element(B, MPO, A) + + +def overlap(A, B): + return mps.overlap(A, B) + +def matrix_element(A, M, B): + return mps.matrix_element(A, M, B) + +def contract(A, B=None, MPO=None): + return mps.contract(A, B=B, MPO=MPO) + diff --git a/engine/heom/mps_utils.py b/engine/heom/mps_utils.py new file mode 100644 index 0000000..e080aca --- /dev/null +++ b/engine/heom/mps_utils.py @@ -0,0 +1,212 @@ +import numpy as np + +def mapint(i, N): + if i < 0: + i += N + if i < 0 or i >= N: + raise IndexError("The index (%d) is out of range." % i) + return i + +def determine_truncation(S, tol = None, nbond = None, chimin = None, is_ortho = True): + nsvd = S.shape[0] + + #we will only allow for truncation if we are currently at the orthogonality centre + if is_ortho: + if not nbond == None: + if nsvd > nbond: + nsvd = nbond + + if not tol == None: + #print(np.amax(S), tol) + #print(S/np.amax(S)) + ntrunc = np.argmax(S < tol*np.amax(S)) + if(ntrunc == 0): + ntrunc = nsvd + if nsvd > ntrunc: + nsvd = ntrunc + + if not chimin == None: + if nsvd < chimin: + nsvd = chimin + + return nsvd + + +def two_site_matrix_to_tensor(M, d1, d2): + return np.transpose(M.reshape((d1, d2, d1, d2)), [0, 2, 1, 3]) + +#functions used for constructing two site object +def two_site_mpo(M, d1, d2, order = 'mpo', tol=None): + Mmt = two_site_matrix_to_tensor(M, d1, d2) + Mm = Mmt.reshape((d1*d1, d2*d2)) + + Q, S, Vh = np.linalg.svd(Mm, full_matrices=False, compute_uv=True) + nsvd = determine_truncation(S, tol=tol) + + Q = Q[:, :nsvd] + S = np.diag(S[:nsvd]) + Vh = Vh[:nsvd, :] + R = Q @ S + if(order == 'mpo'): + Rt = R.reshape((1, d1, d1, nsvd)) + Vht = Vh.reshape((nsvd, d2, d2, 1)) + else: + Rt = np.transpose(R.reshape((1, d1, d1, nsvd)), [0, 1, 3, 2]) + Vht = Vh.reshape((nsvd, d2, 1, d2)) + return Rt, Vht + + +def permute_nsite_dims(M, ind, ds): + indms = [ [x, i] for i,x in enumerate(ind)] + #for i, m in enumerate(mi): + # indms.append([m, i]) + + + #sort the list based on the values of m + indms = sorted(indms, key = lambda x : x[0]) + pind = [ims[0] for ims in indms] + + perms = [ims[1] for ims in indms] + + pdms = [ds[x] for x in perms] + + for ims in indms: + perms.append(ims[1]+len(ds)) + + Mtens = np.transpose(M.reshape((*ds, *ds)), axes=perms) + return Mtens, pind, pdms + + +def nsite_tensor_to_mpo(M, d): + #permute the matrix into a tensor with indices the same as needed for an MPO + xs = [val for pair in zip(range(len(d)), range(len(d), 2*len(d))) for val in pair] + #now we need to permute this into the mpo ordering + return np.transpose(M, axes=xs) + + +def nsite_mpo(M, d, order = 'mpo', tol=None): + Mt = nsite_tensor_to_mpo(M, d) + + mpo = [] + + d1 = 1 + d2 = int(np.prod(d))**2 + + nsvdp = 1 + + #now we need to iterate over each of the indices of Mt and + for i in range(len(d)-1): + d1 = nsvdp * (d[i]**2) + d2 = d2 // (d[i]**2) + Q, S, Vh = np.linalg.svd(Mt.reshape((d1, d2)), full_matrices=False, compute_uv=True) + nsvd = determine_truncation(S, tol=tol) + + Q = Q[:, :nsvd] + S = np.diag(S[:nsvd]) + Vh = Vh[:nsvd, :] + Mt = S @ Vh + + d1 = nsvd + if order == 'mpo': + mpo.append(Q.reshape((nsvdp, d[i], d[i], nsvd))) + else: + mpo.append(np.transpose(Q.reshape((nsvdp, d[i], d[i], nsvd)), [0, 1, 3, 2])) + + if(i + 2 == len(d)): + if(order == 'mpo'): + mpo.append(Mt.reshape(nsvd, d[-1], d[-1], 1)) + else: + mpo.append(Mt.reshape(nsvd, d[-1], 1, d[-1])) + nsvdp=nsvd + return mpo + + +def identity_pad_mpo(nsvd, d, order = 'mpo'): + if(order == 'mpo'): + return np.transpose(np.identity((nsvd*d)).reshape(nsvd, d, nsvd, d), [0, 1, 3, 2]) + else: + return np.identity((nsvd*d)).reshape(nsvd, d, nsvd, d) + + +def check_compatible(A, B, Atype, Btype): + if not (isinstance(A, Atype) and isinstance(B, Btype)): + raise ValueError("Unable to contract two non-mps object.") + if(not A.nsites() == B.nsites()): + raise ValueError("Unable to contract mps object with different number of sites.") + + if(A.local_dimension(slice(0, len(A))) != B.local_dimension(slice(0, len(A)))): + raise ValueError("Unable to contract mps objects with different local hilbert space dimensions.") + + +#functions used to shift the orthogonality centre of an mps or mpo +def update_mps_ortho_centre(il, ir, oc, dir): + if dir == 'right': + if oc == il: + return ir + else: + return oc + elif dir == 'left': + if oc == ir: + return il + return oc + return oc + + +def local_canonical_form(Mi, Mj, dir, il, oc, tol = None, nbond = None, chimin = None): + #if dir == 'right' we need to reshape and svd the left tensor + if dir == 'right': + dims = Mi.shape + A = Mi.reshape((dims[0]*dims[1], dims[2])) + + Q = None; R = None + if tol == None and nbond == None: + Q, R = np.linalg.qr(A, mode='reduced') + + else: + Q, S, Vh = np.linalg.svd(A, full_matrices=False, compute_uv=True) + + nsvd = determine_truncation(S, tol = tol, nbond = nbond, chimin=chimin, is_ortho = (il == oc)) + + + Q = Q[:, :nsvd] + S = np.diag(S[:nsvd]) + Vh = Vh[:nsvd, :] + R = S @ Vh + + return Q.reshape((dims[0], dims[1], R.shape[0])), R, Mj + + #otherwise we reshape and svd the right tensor + elif dir == 'left': + dims = Mj.shape + B = Mj.reshape((dims[0], dims[1]*dims[2])) + + Vh = None; R = None + if tol == None and nbond == None: + U, L = np.linalg.qr(B.T, mode='reduced') + Vh = U.T + R = L.T + + else: + Q, S, Vh = np.linalg.svd(B, full_matrices=False, compute_uv=True) + + nsvd = determine_truncation(S, tol = tol, nbond = nbond, chimin=chimin, is_ortho = (il + 1 == oc)) + + Q = Q[:, :nsvd] + S = np.diag(S[:nsvd]) + Vh = Vh[:nsvd, :] + R = Q @ S + + return Mi, R, Vh.reshape( (R.shape[1], dims[1], dims[2])) + + else: + ValueError("Invalid dir argument") + +def shift_mps_bond(Mi, Mj, dir, il, oc, tol = None, nbond = None, chimin = None): + X, R, Y = local_canonical_form(Mi, Mj, dir, il, oc, tol=tol, nbond=nbond, chimin=chimin) + if dir == 'right': + return X, np.tensordot(R, Y, axes=([1], [0])) + #otherwise we reshape and svd the right tensor + elif dir == 'left': + return np.tensordot(X, R, axes=([2], [0])), Y + else: + ValueError("Invalid dir argument") diff --git a/heom_dynamics.py b/heom_dynamics.py new file mode 100644 index 0000000..d373794 --- /dev/null +++ b/heom_dynamics.py @@ -0,0 +1,138 @@ +from engine.heom import * + +import numpy as np +import scipy as sp +import h5py + +import sys + +#current function used to define the support points for the aaa algorithm. The best choice of support points will depend on the nature of the spectral function you are considering (e.g. where are the interesting features of the spectral function, discontinuities, derivative discontinuities, sharp peaks) +#This choice works rather well for exponential cutoffs with challenging points at zero but won't work well for generic spectral densities. + +def heom_dynamics(rho0, H, baths, dt, tmax, chimax, tol = 1e-6, plot_aaa = False, fname=None, output_stride =10): + nhilb = rho0.shape[0] + nliouv = nhilb*nhilb + + bs = [x.discretise(output_fitting = plot_aaa) for x in baths] + + if plot_aaa: + import matplotlib.pyplot as plt + for b in bs: + plt.plot(b["wf"], b["Sw"]) + plt.plot(b["wf"], b["Sw_fit"]) + plt.show() + + nt = int(tmax/dt)+1 + + #now build the propagator matrices for HEOM + print("Build HEOM propagator matrices", file=sys.stderr) + prop = HEOM_propagator(bs, dt) + + A = setup_HEOM_ados(rho0, prop) + + print("Build HEOM propagator MPO", file=sys.stderr) + + rhos = np.zeros((nt+1, nliouv), dtype = np.complex128) + + rhos[0, :] = rho0.flatten() + + print(0, end=' ' ) + for j in range(rhos.shape[1]): + print(np.real(rhos[0, j]), np.imag(rhos[0, j]), end=' ') + print(A.maximum_bond_dimension()) + + n = 1 + for i in range(nt): + #update the time dependent system propagator at the midpoint of the next step + Hsys = H(i*dt+dt/2) + Lsys = commutator(Hsys) + Us = sp.linalg.expm(-1.0j*dt*Lsys) + + #apply the propagator through a time step + apply_propagator(Us, prop, A, tol=tol, nbond=chimax) + + #compute rho from the object + rho = extract_rho(A) + + #and save it in the rhos array + rhos[i+1, :] = rho + print((i+1)*dt, end=' ' ) + for j in range(len(rho)): + print(np.real(rho[j]), end=' ') + print(np.real(rho[0]+rho[3]), end=' ') + print(A.maximum_bond_dimension()) + sys.stdout.flush() + + if (i % output_stride == 0 and not fname == None): + h5out = h5py.File(fname, 'w') + h5out.create_dataset('t', data=np.arange(nt+1)*dt) + h5out.create_dataset('rho', data=rhos) + h5out.close() + + if not fname == None: + h5out = h5py.File(fname, 'w') + h5out.create_dataset('t', data=np.arange(nt+1)*dt) + h5out.create_dataset('rho', data=rhos) + h5out.close() + return rhos + + +def main(): + sx = np.array([[0, 1], [1, 0]], dtype = np.complex128) + sy = np.array([[0, -1.0j], [1.0j, 0]], dtype=np.complex128) + sz = np.array([[1, 0], [0, -1]], dtype = np.complex128) + + #and the system part of the system bath coupling operator + S = sz + + #function defining the bath spectral function + beta = 5 + s = 1 + zeta = 0.5 + wc = 7.5 + def Jw(w): + return np.sign(w)*np.pi/2.0*zeta*wc*np.power(np.abs(w)/wc, s)*np.exp(-np.abs(w)/wc) + + def Sw(w): + if beta == None: + return Jw(w)*np.where(w > 0, 1.0, 0.0) + else: + return Jw(w)*0.5*(1+1.0/np.tanh(beta*w/2.0)) + + + #the initial value of the system density operator + rho0 = (np.identity(2) + sz)/2 + + #set up the evolution parameters + dt = 0.05 + tmax = 15.0 + + #setup the ados parameters - accuracy of representation of auxiliary density operator + nbose = 15 + Lmin = 5 + tol = 1e-8 + chimax = 100 + + #setup the aaa parameters - controls accuracy of representation of bath correlation function + #convergence with respect to the aaa_tol parameter depends significantly on the choice of support points used, and also on the form of the bath correlation function. + #Here I am using a softmspace grid with a rather tight zero. This approach seems to work quite well for baths with an exponential cutoff at zero temperature but doesn't necessarily work the best for e.g. Debye baths. + #The choice of support points is somewhat an art and the accuracy can really only be determined by monitoring both the accuracy of the aaa fitting of the spectral function and also the correlation function. + #To change the + aaa_tol = 1e-3 + wmax = wc + + #construct the HEOM baths object + baths = [heom_bath(S, Sw, nbose, Lmin=Lmin, wmax = wc, aaa_tol=aaa_tol)] + + eps = 0.0 + delta = 1.0 + + def H(t): + return eps*sz + delta*sx + + + rhos = heom_dynamics(rho0, H, baths, dt, tmax, chimax, tol=tol, plot_aaa=True, fname='sbm_dynamics.h5') + +if __name__ == "__main__": + main() + diff --git a/visualise_heom.py b/visualise_heom.py new file mode 100644 index 0000000..2caeea9 --- /dev/null +++ b/visualise_heom.py @@ -0,0 +1,38 @@ +import numpy as np +import matplotlib.pyplot as plt +import h5py +import argparse + +def plot(fname, inds): + h5 = None + try: + h5 = h5py.File(fname, 'r') + t = np.array(h5.get('t')) + rho = np.array(h5.get('rho')) + + h5.close() + except: + print("Failed to read input file") + return + + + for i in inds: + try: + plt.plot(t, np.real(rho[:, i]), '-', label=r'$\mathrm{Re}(\rho_{%d}(t))$'%(i)) + plt.plot(t, np.imag(rho[:, i]), '--', label=r'$\mathrm{Im}(\rho_{%d}(t))$'%(i)) + except: + print("Failed to plot index: %d"%(i)) + plt.legend(frameon=False) + plt.show() + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description='Plot density matrix elements output by a heom calculation.') + parser.add_argument('fname', type=str) + parser.add_argument('--inds', nargs='+', default = [0]) + + args = parser.parse_args() + x = [int(i) for i in args.inds] + plot(args.fname, x) + + From 99786bedfc193d5f7353227acaca896ac510cf0b Mon Sep 17 00:00:00 2001 From: Lachlan Lindoy Date: Fri, 16 Aug 2024 11:05:28 +0100 Subject: [PATCH 2/2] Added git ignore. --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7bb7933 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.py[cod] +__pycache__/