From ff14e6f94ef1f82cb098c0aec72a1e137ced6577 Mon Sep 17 00:00:00 2001 From: Vasil Saroka Date: Sat, 10 Apr 2021 15:46:37 +0300 Subject: [PATCH] Sneg Sneg library by Rok Zitko has been added. It realizes second quantization technique. --- PacletInfo.m | 6 +- TBpack/CustomTicks.m | 12 +- .../Symbols/ElectronicStructure.nb | 65 +- TBpack/Logo_3.png | Bin 0 -> 4519 bytes TBpack/MaTeX.m | 6 +- TBpack/Sneg.m | 4271 +++++++++++++++++ TBpack/TBpack.m | 6 +- 7 files changed, 4317 insertions(+), 49 deletions(-) create mode 100644 TBpack/Logo_3.png create mode 100644 TBpack/Sneg.m diff --git a/PacletInfo.m b/PacletInfo.m index 7ff18f9..bbc017b 100644 --- a/PacletInfo.m +++ b/PacletInfo.m @@ -1,10 +1,10 @@ (* Paclet Info File *) -(* created 2021/03/04*) +(* created 2021/04/08*) Paclet[ Name -> "TBpack", - Version -> "0.3.0", + Version -> "0.4.0", MathematicaVersion -> "10.0+", Description -> "Tight-binding calculations in Mathematica.", Creator -> "Vasil Saroka <40.ovasil@gmail.com>", @@ -23,7 +23,7 @@ Publisher -> "vasilsaroka", URL -> "https://github.com/vasilsaroka/TBpack", Support -> "https://t.me/vasilsaroka", - Icon -> "Logo_2.png", + Icon -> "Logo_3.png", UUID -> "MaTeX documentation embeded", Location -> "https://github.com/vasilsaroka/TBpack", Extensions -> diff --git a/TBpack/CustomTicks.m b/TBpack/CustomTicks.m index fc57a0b..0b5256e 100644 --- a/TBpack/CustomTicks.m +++ b/TBpack/CustomTicks.m @@ -1,12 +1,16 @@ (* ::Package:: *) (* :Title: CustomTicks *) -(* :Context: CustomTicks` *) (* :Author: Mark A. Caprio, Department of Physics, University of Notre Dame *) +(* :Context: CustomTicks` *) +(* :Version: 2.1.0 *) +(* :Date: 2016-03-12 *) + +(* :Mathematica Version: 6.0+ *) +(* :Copyright: (c) 2016, Mark A. Caprio *) + + (* :Summary: Custom tick mark generation for linear, log, and general nonlinear axes. *) -(* :Copyright: Copyright 2016, Mark A. Caprio *) -(* :Package Version: 2.1.0 *) -(* :Mathematica Version: 6.0 *) (* :History: MCAxes package, January 10, 2003. MCAxes and then MCTicks packages distributed as part of LevelScheme, 2004. diff --git a/TBpack/Documentation/English/ReferencePages/Symbols/ElectronicStructure.nb b/TBpack/Documentation/English/ReferencePages/Symbols/ElectronicStructure.nb index 9dbd313..4f57470 100644 --- a/TBpack/Documentation/English/ReferencePages/Symbols/ElectronicStructure.nb +++ b/TBpack/Documentation/English/ReferencePages/Symbols/ElectronicStructure.nb @@ -10,10 +10,10 @@ NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 158, 7] -NotebookDataLength[ 12110933, 219921] -NotebookOptionsPosition[ 11994842, 217724] -NotebookOutlinePosition[ 11995609, 217750] -CellTagsIndexPosition[ 11995524, 217745] +NotebookDataLength[ 12110600, 219912] +NotebookOptionsPosition[ 11994508, 217715] +NotebookOutlinePosition[ 11995275, 217741] +CellTagsIndexPosition[ 11995190, 217736] WindowFrame->Normal*) (* Beginning of Notebook Content *) @@ -216476,25 +216476,16 @@ Cell[BoxData[ CellLabel->"In[9]:=", CellID->472538745,ExpressionUUID->"8c315538-0ffe-481e-883a-2d4d9272f070"], -Cell[TextData[{ - "Energy levels for several organic molecules which structures are invoked \ -with the built-in function ", - Cell[BoxData[ - ButtonBox["Molecule", - BaseStyle->"Link"]], "InlineFormula",ExpressionUUID-> - "c8e21549-9797-49e2-b2ba-1cf51286106b"], - " (built-in function ", - Cell[BoxData[ - ButtonBox["MoleculePlot", - BaseStyle->"Link"]], "InlineFormula",ExpressionUUID-> - "cac57b12-50a1-4a13-892a-ec9851a41caf"], - " is used for their structure visualization in the inset):" -}], "ExampleText", +Cell["\<\ +Energy levels for several organic molecules and their structure visualization \ +in the inset:\ +\>", "ExampleText", CellChangeTimes->{{3.807296519639682*^9, 3.807296578926073*^9}, { 3.8072966132410355`*^9, 3.807296652813299*^9}, {3.8161814906090684`*^9, 3.8161815604261713`*^9}, {3.816181597990226*^9, 3.816181599741228*^9}, { 3.8161816357252817`*^9, 3.8161816419722905`*^9}, {3.816182088673994*^9, - 3.8161820977590075`*^9}, {3.8161821388460665`*^9, 3.816182176278123*^9}}, + 3.8161820977590075`*^9}, {3.8161821388460665`*^9, 3.816182176278123*^9}, { + 3.8259454422758574`*^9, 3.825945497609022*^9}}, CellID->83758257,ExpressionUUID->"3ab175fe-838a-401f-8d67-d86abc5c72c3"], Cell[CellGroupData[{ @@ -217226,7 +217217,7 @@ AMy+5XY= CellLabel->"Out[10]=", CellID->42744818,ExpressionUUID->"1570f5f5-a5fd-48a7-bfa4-14be890d547c"] }, Open ]] -}, Closed]], +}, Open ]], Cell[CellGroupData[{ @@ -217744,7 +217735,7 @@ CellTagsIndex->{ *) (*CellTagsIndex CellTagsIndex->{ - {"ExtendedExamples", 11995329, 217738} + {"ExtendedExamples", 11994995, 217729} } *) (*NotebookFileOutline @@ -219885,41 +219876,41 @@ Cell[11928861, 216183, 12786, 275, 1446, "Output",ExpressionUUID->"a5d5f469-652c Cell[CellGroupData[{ Cell[11941708, 216465, 327, 6, 26, "ExampleSection",ExpressionUUID->"dc8a62bd-acee-4765-bc19-53baa96a786f", CellID->1653164318], -Cell[11942038, 216473, 170, 3, 37, "Input",ExpressionUUID->"8c315538-0ffe-481e-883a-2d4d9272f070", +Cell[11942038, 216473, 170, 3, 39, "Input",ExpressionUUID->"8c315538-0ffe-481e-883a-2d4d9272f070", CellID->472538745], -Cell[11942211, 216478, 951, 19, 52, "ExampleText",ExpressionUUID->"3ab175fe-838a-401f-8d67-d86abc5c72c3", +Cell[11942211, 216478, 617, 10, 30, "ExampleText",ExpressionUUID->"3ab175fe-838a-401f-8d67-d86abc5c72c3", CellID->83758257], Cell[CellGroupData[{ -Cell[11943187, 216501, 19718, 479, 982, "Input",ExpressionUUID->"8849a9db-de1b-490c-89c1-42ac1773d038", +Cell[11942853, 216492, 19718, 479, 1009, "Input",ExpressionUUID->"8849a9db-de1b-490c-89c1-42ac1773d038", CellID->66424741], -Cell[11962908, 216982, 11259, 244, 428, "Output",ExpressionUUID->"1570f5f5-a5fd-48a7-bfa4-14be890d547c", +Cell[11962574, 216973, 11259, 244, 428, "Output",ExpressionUUID->"1570f5f5-a5fd-48a7-bfa4-14be890d547c", CellID->42744818] }, Open ]] -}, Closed]], +}, Open ]], Cell[CellGroupData[{ -Cell[11974216, 217232, 249, 5, 26, "ExampleSection",ExpressionUUID->"98823c6f-0659-41a6-a97c-3dd905663f98", +Cell[11973882, 217223, 249, 5, 46, "ExampleSection",ExpressionUUID->"98823c6f-0659-41a6-a97c-3dd905663f98", CellID->589267740], -Cell[11974468, 217239, 335, 6, 45, "ExampleText",ExpressionUUID->"fc7b6cff-b06e-4d72-a2cf-5ac820df27a1", +Cell[11974134, 217230, 335, 6, 45, "ExampleText",ExpressionUUID->"fc7b6cff-b06e-4d72-a2cf-5ac820df27a1", CellID->158555077], -Cell[11974806, 217247, 169, 3, 55, "Input",ExpressionUUID->"6d1e7201-b832-456a-8fe4-d84b45a81fdb", +Cell[11974472, 217238, 169, 3, 55, "Input",ExpressionUUID->"6d1e7201-b832-456a-8fe4-d84b45a81fdb", CellID->19388544], -Cell[11974978, 217252, 288, 5, 45, "ExampleText",ExpressionUUID->"15d5ab26-c0c4-4483-b137-18af9b8251e3", +Cell[11974644, 217243, 288, 5, 45, "ExampleText",ExpressionUUID->"15d5ab26-c0c4-4483-b137-18af9b8251e3", CellID->225468900], -Cell[11975269, 217259, 2117, 53, 381, "Input",ExpressionUUID->"e89a2b79-ac05-4999-bfa8-75c375a2f297", +Cell[11974935, 217250, 2117, 53, 381, "Input",ExpressionUUID->"e89a2b79-ac05-4999-bfa8-75c375a2f297", CellID->100735450], -Cell[11977389, 217314, 418, 7, 75, "ExampleText",ExpressionUUID->"f8afc424-a463-45ab-8782-b24ac9b5014b", +Cell[11977055, 217305, 418, 7, 75, "ExampleText",ExpressionUUID->"f8afc424-a463-45ab-8782-b24ac9b5014b", CellID->230233141], -Cell[11977810, 217323, 1669, 38, 178, "Input",ExpressionUUID->"0effa785-0442-42cf-969e-5b939e650d6f", +Cell[11977476, 217314, 1669, 38, 178, "Input",ExpressionUUID->"0effa785-0442-42cf-969e-5b939e650d6f", CellID->774095697], -Cell[11979482, 217363, 765, 17, 109, "ExampleText",ExpressionUUID->"9f8e36ec-807c-44c8-86cb-72b03687f065", +Cell[11979148, 217354, 765, 17, 109, "ExampleText",ExpressionUUID->"9f8e36ec-807c-44c8-86cb-72b03687f065", CellID->67501718], Cell[CellGroupData[{ -Cell[11980272, 217384, 5226, 132, 705, "Input",ExpressionUUID->"d04b97fc-928f-46bf-a350-a2397dedb4f2", +Cell[11979938, 217375, 5226, 132, 705, "Input",ExpressionUUID->"d04b97fc-928f-46bf-a350-a2397dedb4f2", CellID->38364410], -Cell[11985501, 217518, 7499, 160, 835, "Output",ExpressionUUID->"5c8da57b-5af4-4fd9-abfd-1fbbef8b41d8", +Cell[11985167, 217509, 7499, 160, 835, "Output",ExpressionUUID->"5c8da57b-5af4-4fd9-abfd-1fbbef8b41d8", CellID->32333434] }, Open ]], -Cell[11993015, 217681, 1799, 39, 174, "ExampleText",ExpressionUUID->"29e836c5-9888-49ea-82bb-44a1273c1b93", +Cell[11992681, 217672, 1799, 39, 174, "ExampleText",ExpressionUUID->"29e836c5-9888-49ea-82bb-44a1273c1b93", CellID->45530709] }, Closed]] }, Open ]] diff --git a/TBpack/Logo_3.png b/TBpack/Logo_3.png new file mode 100644 index 0000000000000000000000000000000000000000..3a4c0d2cf93d7163ffad8fea34572684000e9729 GIT binary patch literal 4519 zcmXAt3p`W*|Hp^L+_JK?G)yh6a&InaVKPi^6EabTVeVuhxeY6GsraIX#i|veEcA7W z$}KI?cgZc}J2jDe?y3Lz{{H9jIOjYb=W#Bd^LfACuh;XO3?~O`c^PFH5C|ksu(5Oo z_PYPA5OH8T=pEJv>>gfmagTN;#~>r3!u_cs6l8Q<1O-WnrTT+FvBR#8Bl7J?xvc=^ zesvuLT&m8cQ9kU~@lcwA>`m%mH%j{PfVR-{TwHlq%S{u1UHy_V|wM= z7138^b)S0%q_V_ozYbvK_pop~^X#HGt*ouIGz=oWCEOnk+WXvd@Z-s61pVu4pgs8R z&$P5CL1Y#7-oPutN$A%S#ZCSErlsZV80Nzf}yxDux3G zvC(+<=(FMe(J|yG3g~ncg&a-sL&i{}1Ch~z6y&+^5dSE0DAJA`7C<8hP>?3bn3x#t znEyYl?HeBI$>F90_fh}vKL5L%8Xkspp@vdG`nvjhCb}3CU1Oxau8F>}$#H#T@OUr? z1l1#0;#}O*OWbF}UHTixgIinYPwA96oFAs1NpUtW1rw3GEiq{6T`11Pq_Qx;EU3D; z=yq}Nx&)>ZER z?Af!Bva$ayW@Zdk7}WCk@#D+yPtHdm0!A(e1c;vXnl^=vHvX};Bg36x`TVbANq@Fi zW(K5iP{zdAi}-=Zq|;Ea#iQ*>iGg(UUiPW2e=~JjR)>SHDX3lcT6=fr__=CNb%H*G z{wEldnVG4U^s5HV%(n8SP&$^eI4IO2#(Mc>vTi4e;9zg>dNepiS3$+3wz>K0?_V47V~t7Xpj`vBwfS*1JQ2>Qp7}Jhaw>jxs*!D{q2tvU zy)@Z9@;ux~S$=?)icH4Ocbj#t|NJhPf3Y@p@CKF_JWy?UU>BZP+0W$I6nvGC+WBjV zud@cW95+Iy-4q&Y@;iNMbmjIjr0K{K2j8q4h!at(KZYhICVXmQfdN#Vf6gF3q7HAa z%t+JSJw5&4+@TPjW+ipFy}ey@H!1Pw)CaBe^{6}oi+dvHlGnKc3Jx0WHGr9A`+h!* zoYE0sTwT*c_Aw%7hwFgJeBgW}N4QfJyS_RS9z&xY+MVI2f5Obn%+z$Qv3h=s&*ztv zl+=eF-klMPH8nI`ZC7;sG1_EQ;?VYQf6&UOs`H4h(D9a2a&7VYSN`Oxr|EW*;g@N& zSHN>qk%5LqU|bscwhbyr<)lB$->|4hw9MpWWo6w&4TTu)&$fay;6gT=EqLf?q?}!N zF_9#wMKg6fi7WNw&l z2ErLA!cnFL)KVjlR2S6Wk-0ZMK3=8kM9duG{=~z(iMc51^0yw0y&_nh`R|)IZ&(Jw zx}AN03MhnJaFMC*X);_}N9TF`$XbE`Gg0SUjsf01j*TLurlxlL{X_E_pdHJ%ZSYiu z!{3OOb5XNnFSe#W+~c5G+&jIZXG0-SFhm_s^r$MJkRWB{X}QN}zfo9oz)buZ$?~G&Vn31_5w_k~TWgs+nSbn@ zm6a9a+O<=?iuhZo?pLo~T{WHGN<5m+p*A)(#V01VJbY*k!!tQSb;=Pn6%{!`Oi@oz z9dLE~a`XHagN!h}YMgHiw3?s(#>-@=1rgTs{{8Ces>W_Wi3Np)U7ekhl9HdQz4*)6 zDXTh<*_oO1o}OvLdUC0Hfb%@n>%2`Nxw;-<%CtSb@MIZ#*Y?5E#?^ZtAaw7rFm0vG z+0PBj)BTkUa=VY-7g$=~_VN4Zx7KP2wd6L2%TWV8y_&e?ZgczR4ISpy3K1Ohc zGl;ON3olrx?x==hhlimfVHoA)O`sz|w$|2Cs0`%wmX;PkJ+=?TIVh%&{+(Z&Gf(8( zVwbuP`_>F5h>KeuWAbP-pMcWz_)|y)lzr`BUr$XXpq&5PV0(>;}3!f=X2SpU$_zRV6Tm+SBE?@(e=cX@W$tLbui#b%d8c+1AO zQL%J#AW4bOqa$INcv69IasBOEW<^CVsz&oX8s=Q?qNx;YShR7~RPHnm`dT*{VLbq7 zDy0pK8?2?C)l^ahi}j7JLU2$pu5_l649~ayy1po+BW-fQVxu$HIOqzu%9SfGhC(VT z043>l@-!tfEs1_CKvW#`Zm`~|Q>Ts}SIq3l?e=@%d#uM{bk&_( z#r8%b7-Tp@&;BeFEWjj;i;QzbaX9E3vn|c8NSoZOE+Jb)gdKv2+Y7i<7FYTq48OQG zo+M>8H9xPWqLPj_us`qNp)V(4Qt5J(qQ5UkD3flXT?mx_aeojDABpJVY2u)AZLw_8 zmT9rP-cY_Qu-eq6oRdzrTO~?y23SD`x?@avloK5pvLP z{q*ZYm$6S~u2l9= z$Y^SJgmt)88elMh#te%RwpRvpI}aW_xG-)}=FRf2^DwXnOyDS!LmfVR_v+H3+3yMC z9LmV(=qR#=EgGtq|Cyr47|`N3MtzmR6It9_G0SCSxRW=Fy94Ls>MFpzd@MBHxH=YD z9Q@~(D5LfF%0}z2k-_!#@YZcf9R@jY#1K#H_TxFBonT+%=90Fz=i?Qyc^9jTjUzcO zM}b)6@eVM3Vq%_CS!@yR3P^6GJ`|Ypl`Gp0-GD;yL_zKQ_YS?wbYv-7#``@l1kUK| zf5$SA*xokk3h$rJ?ls*C5peAnd^GYF695MSGDc%}M&@3f-%+z30pV)deM!G(!+})m z^?%0VlCw_aEGOLMQlF-n#}i>?GaKApFl z^6jM|M0}lA7+$Z$p{e>9!~(18_)q}DXJ=>64m~?8xU+d9DQUarRJ?ycfLd_2Plp~PF$amA{7{gkIC^J936|eCyb2)Np2=3JLxZZn%6tm z64IcA+$y(YN||N#!WvD$-LfKITrthY*WWx9px5^AlI+r@WhBMx~ zckdy>Zg80sI+{k?fkQ&*NJxEX9!i>qgY?g=IG2AMA0H1foM_FqlCYE>wlAQVLkLjI z2hJo<&-+y6VFT40d>+=!ENN>sbb!{>)Fe$u=i9>YN0{Ot`$7zlFnvgddEIEHcDHufDYtkLB)Lm zV3R((=ekzsAZ--Q+}zx>Bm$3*h5VBE^@^uy0dfC8^Lmp)Wx?Nq0KUk##S3cPDmx+(_JXGcGmRZSH9!%nVw@n?;nT8&CG3;R zO4$hcsl@~}B`QYw<~kDM9U;KTTV!zFL?Q-gRF<3T9vN#WSmj1n%I_mU=Qp}lesA3Q z^&`!zfYKESvMU`d)x!3s&Mjfb2L_DH%r0sd9wGqTncd;DHYhc9kja5F#CSW}lfx(C z*g|7tkeJ5B0?NYc(#$@~#%gOKSbQ&@sHr4jC`&gNU?Px^GQbR00fYwm>M#%Rc_%Ob zh@EQ?#j>R)7qO0l+L^|_ji-7@09JBQDqUiC!D6>SJ6v6Htx$_(h4HnNp;DQZE2ELc z)z#H#rncwb-f&3^bWso6gvR842NEYaI%& zfKZ^y$q{uP`ILRSyYL?>Tr0!Inm5*kH`+eCm%9KX*FABf;px+bzS9bEY$JVr@r+g= zk5^|$V&=!>TNK(s+^C4o8uxc=+deh!PX?qUCFf#;z5q7$D(duuSDs3l_$h+Dz{M3P z1P}vsq~9|V4$4mY4aA|w;llug#U&;tvdOmt*m06@KQbAB6Zy~uOe|AUg z0S^b&{qw!HqQd)=HwkK+Ba~}P=0oA@WP^xwATL$=&iZ1)o#3y?(_HcbX@h?<&WHInz0mP$s!Pe(Q^AvC2 z0*aXX!pUEQ8nV|VV3OdSwc_s8qa+6TWx8&b91&JtA(JNp0FTNlcahQ7omX^<-;WF{ zy*#z^=HZcgoRl|0auApI`WhMTZIA|)3$716E+>JqhFW4(uh~hE6y)XW@>hO~)wPND zy_{ZfXc{mV16(>nuQNyZGM_Z!_PpE5wLhxiU1Tu_ZBE75x79rZ(4}HEeI-dkgh}6% UhX4fue;Pmpyn|)6g*W~G0D+oS;{X5v literal 0 HcmV?d00001 diff --git a/TBpack/MaTeX.m b/TBpack/MaTeX.m index 2d6ac7a..ba8c7a6 100644 --- a/TBpack/MaTeX.m +++ b/TBpack/MaTeX.m @@ -1,13 +1,15 @@ (* Mathematica Package *) (* Created by IntelliJ IDEA http://wlplugin.halirutan.de/ *) +(* ::Package:: *) + (* :Title: MaTeX *) (* :Author: Szabolcs Horvat *) (* :Context: MaTeX` *) (* :Version: 1.7.8 *) (* :Date: 2015-03-04 *) -(* :Mathematica Version: %%mathversion%% *) +(* :Mathematica Version: 10.0+ *) (* :Copyright: (c) 2019 Szabolcs Horvat *) @@ -246,7 +248,7 @@ $configOK = pdflatexOK && gsOK && cacheSizeOK && workingDirectoryOK; If[Not[$configOK] && $Notebooks, - Print@StringForm["`` for documentation on configuring MaTeX.", Hyperlink["Click here", "paclet:MaTeX/tutorial/ConfiguringMaTeX"]] + Print@StringForm["`` for documentation on configuring MaTeX.", Hyperlink["Click here", "paclet:TBpack/tutorial/ConfiguringMaTeX"]] ]; ] diff --git a/TBpack/Sneg.m b/TBpack/Sneg.m new file mode 100644 index 0000000..309b481 --- /dev/null +++ b/TBpack/Sneg.m @@ -0,0 +1,4271 @@ +(* ::Package:: *) + +(* :Title: Sneg *) +(* :Author: Rok Zitko *) +(* :Context: Sneg` *) +(* :Version: 1.250.0 *) +(* :Date: 2016-10-15 *) + +(* :Mathematica Version: 6.0+ *) +(* :Copyright: (c) 2016 Rok Zitko *) + +(* + SNEG - Mathematica package for calculations with non-commuting + operators of the second quantization algebra + + Copyright (C) 2002, 2006-2016 Rok Zitko + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + + Contact information: + Rok Zitko + F1 - Theoretical physics + "Jozef Stefan" Institute + Jamova 39 + SI-1000 Ljubljana + Slovenia + + rok.zitko@ijs.si +*) + +(* +A few NAMING CONVENTIONS used in this code for function arguments and patterns: + + z - general complex constant/variable or general Grassman constant/variable + x - general real constant/variable + op - operators + a,b,c - general expressions + i - first index == type, AN=1 for annihilation, CR=0 for creation + j - other indexes, site, spin etc. + k - operator indexes, including the first (type) index + - an exception are SEA ordered operators, where k is used to + denote the second (wavenumber) index [for convenience] + sp, sigma - spin, UP=1 or DO=0 + v - state in number representation + + l - list + i - integer + +Every Fermion operator has at least one index, which specifies whether it +is a creation or annihilation operator. Therefore a generic pattern +is op_[i_, j___]. Moreover, if the operator contains a spin index, +it is by convention always the last index! +*) + + +BeginPackage["TBpack`Sneg`"]; + +(* Usage messages *) + +(* for Functions *) +block::usage = "block[expr] is a wrapper for expressions to prevent +automatic expansion. By definition, nc[block[expr]]=expr. Note also +that the expression withing the block will be evaluated in the usual +way. Use block, for example, to prevent automatic distribution over +Plus."; + +(* Replace indexes with abstract indexes ndxfunc[i] with i=1,...,nrindexes *) +(* Be careful: this is only useful for fine-tuned applications to address +very specific problems! *) +sumAbstractIndex::usage = +"sumAbstractIndex[expr] replaces sum indexes with abstract indexes. +This is useful for simplifications in the case of several equivalent +summation indexes."; +sumNameIndex::usage = +"sumNameIndex[expr] is to be used in conjunction with +sumAbstractIndex[expr]."; + +(* Spin X, Y and Z operators for an arbitrary generalized creation operator +g, i.e. some interesting combination of one-site creation operators. g must +be a pure function of spin index. *) +spinxyzgen::usage = +"spinxyzgen[g] generates x, y, and z components of a S=1/2 spin +operator corresponding to an (operator) expression g. g must be a +pure function of the spin index. For example, the conventional spin +operators are obtained by spinxyzgen[op[CR, #]&]."; + +(* Interleave lists *) +interleave::usage = "interleave[{a1, a2, ...}, {b1, b2, ..}] returns {a1, b1, a2, b2,...}."; + +(* Get sublists with a given key *) +keyselect::usage = "keyselect[list, key] selects sublists with the first element equal to key."; + +(* Merge sublists with the same first element (key) *) +(* Below there is another implementation mergebasis[] which does the same! *) +keymerge::usage = "keymerge[list1, list2] merges two lists, so that the sublists with +the same first element are combined."; + +orthogbasisop::usage = "orthogbasisop[bz1, bz0] orthogonalizes the basis bz1 with respect +to the original basis bz. Both basis sets must be given in the creation +operator representation."; + +orthogbasisvc::usage = "orthogbasisvc[bz1, bz0] orthogonalizes the basis bz1 with respect +to the original basis bz. Both basis sets must be given in the occupation +number representation."; + + + +(* Options *) +PrettyOutput::usage = +"Setting PrettyOutput to False prior to loading the package sneg disables +the pretty printing of sneg expressions."; +PrettyInput::usage = +"Setting PrettyInput to False prior to loading the package sneg disables +the translation of certain markup into sneg functions."; + + + +(* Constants *) +$SnegVersion::usage = "$SnegVersion is the version of sneg package used."; +CR::usage = "CR is a Sneg constant equal to 0; c[CR] is a creation operator."; +AN::usage = "AN is a Sneg constant equal to 1; c[AN] is an annihilation operator."; +DO::usage = "DO is a Sneg constant equal to 0; c[CR, DO] is a spin-down (creation) operator."; +UP::usage = "UP is a Sneg constant equal to 1; c[CR, UP] is a spin-up (creation) operator."; + +PROJ0::usage = "PROJ0 is a Sneg constant equal to 1. It is an argument for projector[] specifying projection to a subspace with zero occupancy."; +PROJUP::usage = "PROJUP is a Sneg constant equal to 2. It is an argument for projector[] specifying projection to a subspace with spin up particle."; +PROJDO::usage = "PROJDO is a Sneg constant equal to 3. It is an argument for projector[] specifying projection to a subspace with spin down particle."; +PROJ2::usage = "PROJ2 is a Sneg constant equal to 4. It is an argument for projector[] specifying projection to a subspace with double occupancy."; +PROJ1::usage = "PROJ1 is a Sneg constant equal to 5. It is an argument for projector[] specifying projection to a subspace with single occupancy."; +PROJ02::usage = "PROJ02 is a Sneg constant equal to 6. It is an argument for projector[] specifying projection to a subspace with either zero or double occupancy."; + + + + +(**** Usage messages ****) +(* Add a link to the relevant documentation page. *) +SetAttributes[UsageWithMore, HoldFirst]; +UsageWithMore[symbol_, string_, tag_:Null] := (symbol::usage = string); + +UsageWithMore[snegfermionoperators, +"snegfermionoperators[a, b, c] declares a, b, and c to be +fermionic operators. They obey canonical anti-commutation relations."]; +UsageWithMore[snegbosonoperators, +"snegbosonoperators[a, b, c] declares a, b, and c to be +bosonic operators. They obey canonical commutation relations."]; +UsageWithMore[snegmajoranaoperators, +"snegmajoranaoperators[a, b, c] declares a, b, and c to be +real (Majorana) fermionic operators."]; +UsageWithMore[snegrealconstants, +"snegrealconstants[x, y] declares x and y to be real constants. +Constants are factored out from operator strings within nc[] blocks. +They are invariant under the conjugation."]; +UsageWithMore[snegcomplexconstants, +"snegcomplexconstants[z, w] declares z and w to be complex constants. +Constants are factored out from operator strings within nc[] blocks."]; +UsageWithMore[sneggrassmanconstants, +"sneggrassmanconstants[z, w] declares z and w to be anti-commuting +Grassman quantities."]; +UsageWithMore[snegfreeindexes, +"snegfreeindexes[k, sigma] declares k and sigma to be indexes. This is +taken into account when simplifying sum[] blocks."]; +UsageWithMore[nc, +"nc[...] represents non-commutative multiplication."]; +UsageWithMore[fermionQ, +"fermionQ[c] returns True if c is a fermionic operator."]; +UsageWithMore[bosonQ, +"bosonQ[c] returns True if c is a bosonic operator."]; +UsageWithMore[operatorQ, +"operatorQ[c] returns True if c is an operator."]; +UsageWithMore[ordering, +"ordering[c] = EMPTY | SEA determined how the fermionic operators +are reordered by defining the vacuum state that corresponds to +the operator c. EMPTY corresponds to an empty state vacuum, while +SEA denotes a Fermi sea filled up to a Fermi level."]; +UsageWithMore[snegOrderedQ, +"snegOrderedQ[op1[i1], op2[i2]] returns True if operators op1[i1] and +op2[i2] are canonically ordered."]; +UsageWithMore[acmt, +"acmt[op1[i1], op2[i2]] defines the value of the anti-commutator between +op1[i1] and op2[i2]. By default, the value is given by the canonical +anti-commutator for fermionic operators."]; +UsageWithMore[cmt, +"cmt[op1[i1], op2[i2]] defines the value of the commutator between +op1[i1] and op2[i2]. By default, the value is given by the canonical +commutator for bosonic operators."]; +UsageWithMore[grassmanQ, +"grassmanQ[z] returns True if z is a Grassman constant or variable."]; +UsageWithMore[komutator, +"komutator[a, b] = nc[a, b] - nc[b,a]"]; +UsageWithMore[commutator, +"commutator[a, b] = nc[a, b] - nc[b,a]", "komutator"]; +UsageWithMore[antikomutator, +"antikomutator[a, b] = nc[a, b] + nc[b, a]"]; +UsageWithMore[anticommutator, +"anticommutator[a, b] = nc[a, b] + nc[b, a]", "antikomutator"]; +UsageWithMore[superkomutator, +"superkomutator[a, b, n] = [a, [a, ... [a, b]]], this is order-n nested +commutator."]; +UsageWithMore[supercommutator, +"supercommutator[a, b, n] = [a, [a, ... [a, b]]], this is order-n nested +commutator.", "superkomutator"]; +UsageWithMore[superantikomutator, +"superantikomutator[a, b, n] = {a, {a, ... {a, b}}}, this is order-n nested +anti-commutator.", "superkomutator"]; +UsageWithMore[superanticommutator, +"superanticommutator[a, b, n] = {a, {a, ... {a, b}}}, this is order-n nested +anti-commutator.", "superkomutator"]; +UsageWithMore[dd, +"dd[expr] denotes the double dots surrounding an expresions, i.e. +a normal ordered string"]; +UsageWithMore[normalorder, +"normalorder[expr] normal orders an expression by subtracting its vacuum +expectation value (vev)."]; +UsageWithMore[normalorderwick, +"normalorderwick[expr] normal orders an expression by subtracting its +vacuum expectation value (vev). Wick's theorem is used to compute the vev.", +"normalorder"]; +UsageWithMore[conj, +"conj[expr] conjugates an expression."]; +UsageWithMore[sum, +"sum[expr, {indexes}] denotes a symbolic sum of an expression +over indexes."]; +UsageWithMore[sumThread, +"sumThread threads sum[{expr1, expr2}, {indexes}] over the list in the first +argument."]; +UsageWithMore[sumExpand, +"sumExpand[expr] expands sums over the addends."]; +UsageWithMore[sumCollect, +"sumCollect[expr] collects sums over the addends."]; +UsageWithMore[sumSimplifyKD, +"sumSimplifyKD[expr] simplifies sums. In particular, it correctly +handles KroneckerDelta functions."]; +UsageWithMore[sumFullSimplify, +"sumFullSimplify[expr] attempts to simplify an expression involving +symbolic sums."]; +UsageWithMore[vev, +"vev[expr] returns the vacuum expectation value of an operator expression."]; +UsageWithMore[inner, +"inner[a, b] returns the inner product between a and b using +nc[] multiplication."]; +UsageWithMore[outer, +"outer[a, b] returns the outer (direct) product between a and b +using nc[] multiplication."]; +UsageWithMore[VMV, +"VMV[v1, m, v2] computes the vector-matrix-vector product using nc[] inner +product multiplication."]; +UsageWithMore[number, +"number[c[i]] returns the number operator corresponding to the operator +c[i]. number[c[i], sigma] with sigma=UP | DO returns the number operator +for a given spin projection only."]; +UsageWithMore[hubbard, +"hubbard[c[i]] returns the local (Hubbard) electron-electron repulsion +operator corresponding to operator c[i]."]; +UsageWithMore[hamiltonian, +"hamiltonian[type, ...] is a stub: it returns the Hamiltonian 'type'."]; +UsageWithMore[chargecharge, +"chargecharge[c[i], d[j]] returns the charge-charge repulsion operator +corresponding to operators (sites) c[i] and d[j]."]; +UsageWithMore[PauliX, "Pauli's X matrix"]; +UsageWithMore[PauliY, "Pauli's Y matrix", "PauliX"]; +UsageWithMore[PauliZ, "Pauli's Z matrix", "PauliX"]; +UsageWithMore[PauliPlus, "Pauli's + (spin raising) matrix", "PauliX"]; +UsageWithMore[PauliMinus, "Pauli's - (spin lowering) matrix", "PauliX"]; +UsageWithMore[spinxyz, +"spinxyz[op[j]] returns the x, y, and z components of the S=1/2 spin +operator corresponding to a fermionic S=1/2 operator op[j]."]; +UsageWithMore[spinall, +"spinall[op[j]] returns the x, y, and z components of a S=1/2 spin +operator, as well as the S^2 operator, corresponding to a fermionic +operator op[j].", "spinxyz"]; +UsageWithMore[spinss, +"spinss[op[j]] returns the S^2 total S=1/2 spin operator squared that +corresponds to a fermionic operator op[j].", "spinxyz"]; +UsageWithMore[spinx, +"spinx[op[j]] returns the x component of a S=1/2 spin operator that +corresponds to a fermionic operator op[j].", "spinxyz"]; +UsageWithMore[spiny, +"spiny[op[j]] returns the y component of a S=1/2 spin operator that +corresponds to a fermionic operator op[j].", "spinxyz"]; +UsageWithMore[spinz, +"spinz[op[j]] returns the z component of a S=1/2 spin operator that +corresponds to a fermionic operator op[j].", "spinxyz"]; +UsageWithMore[spinplus, +"spinplus[op[j]] returns the S=1/2 spin raising operator that +corresponds to a fermionic operator op[j].", "spinxyz"]; +UsageWithMore[spinminus, +"spinplus[op[j]] returns the S=1/2 spin lowering operator that +corresponds to a fermionic operator op[j].", "spinxyz"]; +UsageWithMore[spinspin, +"spinspin[a[i], b[j]] returns the spin-spin scalar product between +two S=1/2 spin operators corresponding to fermionic operators a[i] +and b[j]."]; +UsageWithMore[spinspinpm, +"spinspinpm[a[i], b[j]] returns the +- part of the scalar product +between two S=1/2 spin operators corresponding to fermionic operators +a[i] and b[j].", "spinspin"]; +UsageWithMore[spinspinmp, +"spinspinmp[a[i], b[j]] returns the -+ part of the scalar product +between two S=1/2 spin operators corresponding to fermionic operators +a[i] and b[j].", "spinspin"]; +UsageWithMore[spinspinxy, +"spinspinxy[a[i], b[j]] returns the transverse part of the scalar product +between two S=1/2 spin operators corresponding to fermionic operators +a[i] and b[j].", "spinspin"]; +UsageWithMore[spinspinz, +"spinspinxy[a[i], b[j]] returns the longitudinal part of the scalar product +between two S=1/2 spin operators corresponding to fermionic operators +a[i] and b[j].", "spinspin"]; +UsageWithMore[manyspin, +"manyspin[{a[i], b[j],...}] returns the x, y, and z components of the +total spin, as well as the total spin operator squared, corresponding +to a list of operators given as the argument."]; +UsageWithMore[halfintegerQ, +"halfintegerQ[z] returns True if z is an integer or a half-integer."]; +UsageWithMore[spinmatrixX, +"spinmatrixX[S] returns the matrix of the X component of the spin +operator for arbitrary spin S."]; +UsageWithMore[spinmatrixY, +"spinmatrixY[S] returns the matrix of the Y component of the spin +operator for arbitrary spin S.", "spinmatrixX"]; +UsageWithMore[spinmatrixZ, +"spinmatrixZ[S] returns the matrix of the Z component of the spin +operator for arbitrary spin S.", "spinmatrixX"]; +UsageWithMore[spinmatrixP, +"spinmatrixP[S] returns the matrix of the spin-raising operator +for arbitrary spin S.", "spinmatrixX"]; +UsageWithMore[spinmatrixM, +"spinmatrixM[S] returns the matrix of the spin-lowering operator +for arbitrary spin S.", "spinmatrixX"]; +UsageWithMore[direct, +"direct[b1, b2] returns a direct product of two sets of operator +expressions (nc[] multiplication is used)."]; +UsageWithMore[nambu, +"nambu[c[i], nr] returns the Nambu-spinor components corresponding +to a fermionic operator c[i]. The phase is (-1)^nr. The default +value of nr is 0."]; +UsageWithMore[isospinxyz, +"isospinxyz[c[i], nr] returns the isospin operator corresponding +to the fermionic operator c[i]. The phase is (-1)^nr. The default +value of nr is 0."]; +UsageWithMore[isospin, +"isospin[c[i], nr] returns the isospin operator corresponding +to the fermionic operator c[i]. The phase is (-1)^nr. The default +value of nr is 0.", "isospinxyz"]; +UsageWithMore[nnop, +"nnop[op] is used to define the lattice index that determines +the bipartite lattice used for the isospin operators.", "isospinxyz"]; +UsageWithMore[isospinx, +"isospinx[c[i], nr] returns the x component of the isospin operator +corresponding to the fermionic operator c[i]. The phase is (-1)^nr. +The default value of nr is 0.", "isospinxyz"]; +UsageWithMore[isospiny, +"isospiny[c[i], nr] returns the y component of the isospin operator +corresponding to the fermionic operator c[i]. The phase is (-1)^nr. +The default value of nr is 0.", "isospinxyz"]; +UsageWithMore[isospinz, +"isospinz[c[i], nr] returns the z component of the isospin operator +corresponding to the fermionic operator c[i]. The phase is (-1)^nr. +The default value of nr is 0.", "isospinxyz"]; +UsageWithMore[isospinplus, +"isospinplus[c[i], nr] returns the isospin-raising component of the +isospin operator corresponding to the fermionic operator c[i]. +The phase is (-1)^nr. The default value of nr is 0.", "isospinxyz"]; +UsageWithMore[isospinminus, +"isospinminus[c[i], nr] returns the isospin-lowering component of the +isospin operator corresponding to the fermionic operator c[i]. +The phase is (-1)^nr. The default value of nr is 0.", "isospinxyz"]; +UsageWithMore[manyisospin, +"manyisospin[{ops}, {phases}] returns the X, Y, and Z components +of the total isospin operator for a list of fermionic operators. +The phases of operators must be given as the second argument. +The total isospin squared operator is also returned."]; +UsageWithMore[basis, +"basis[op[i]] returns the basis states (creation operators) for +a single site, i.e. for operator op[i]."]; +UsageWithMore[pow, +"pow[expr, i] raises an expression to the i-th power using the non-commutative +multiplication nc."]; +UsageWithMore[hop, +"hop[a[i], b[j]] returns the electron hopping operator between +sites a[i] and b[j]. hop[a[i], b[j], sigma] with sigma=UP | DO +does the same for a single spin projection sigma."]; +UsageWithMore[holehop, +"holehop[h[i], p[j]] returns the electron hopping operator for +the case of a hole operator h[i] and particle operator p[j]."]; +UsageWithMore[spinfliphop, +"spinfliphop[a[i], b[j]] returns the electron hopping with spin flip +operator between sites a[i] and b[j]."]; +UsageWithMore[twohop, +"twohop[a[i], b[j]] returns the two-particle hopping operator +between sites a[i] and b[j]."]; +UsageWithMore[invertspin, +"invertspin[expr] inverts the spin of the operator expression expr. +The spin index must be the last index of each operator appearing +in the expression."]; +UsageWithMore[makebasis, +"makebasis[ops] constructs a basis of single-particle creation operators +for enlisted operators ops."]; +UsageWithMore[goodopQ, +"goodopQ[c[i]] returns True if the operator c[i] appears in the basis +declared by makebasis."]; +UsageWithMore[op2ndx, +"op2ndx[c[i]] returns the index of operator c[i] in the basis declared +by makebasis."]; +UsageWithMore[ndx2op, +"ndx2op[ndx] returns the creation operator corresponding to +the index ndx in the basis declared by makebasis."]; +UsageWithMore[vacuum, +"vacuum[] returns the vacuum vector in the occupation number (vc) +representation for a basis that has been previously defined using makebasis.", +"vacuumlc"]; +UsageWithMore[ap, +"ap[op, vc] applies an operator op to a vector vc in the occupation +number representation."]; +UsageWithMore[vc, +"vc[...] is a wrapper (vector) that holds the information about the occupancy +of the individual sites in the occupation number representation."]; +UsageWithMore[vc2ops, +"vc2ops[expr] transforms vectors in the occupation-number representation +in expression 'expr' to the corresponding operator strings."]; +UsageWithMore[scalarproductvc, +"scalarproductvc[a, b] computes the scalar product between two expressions +consisting of vectors in the occupation number representation."]; +UsageWithMore[scalarproductop, +"scalarproductop[a, b] computes the vacuum expectation value of the +product between operator expressions a and b. The expression on the +left, a, is automatically conjugated."]; +UsageWithMore[braketvc, +"braketvc[a, op, b] computes the braket , where op is an operator +expression, while a and b are some states (in the occupation number +representation)."]; +UsageWithMore[braketop, +"braketop[a, op, b] computes the braket , where op is an operator +expression, while a and b are some states (in the creation operator +representation).", "braketvc"]; +UsageWithMore[expvop, +"expvop[op, a] computes the braket , where op is an +operator expression and a is some state (in the creation operator +representation).", "braketvc"]; +UsageWithMore[expvvc, +"expvvc[op, a] computes the braket , where op is an operator +expression and a is some state (in the occupation number representation).", +"braketvc"]; +UsageWithMore[matrixrepresentationvc, +"matrixrepresentationvc[op, basis] computes the matrix representation +of the operator op in a given basis (occupation number representation). +matrixrepresentationvc[op, basisL, basisR] returns the matrix +representation of an operator in the case of different basis for +left-multiplication and right-multiplication."]; +UsageWithMore[matrixrepresentationvcfast, +"matrixrepresentationvcfast[op, basis] computes the matrix representation +of the operator op in a given basis (occupation number representation).", +"matrixrepresentationvc"]; +UsageWithMore[matrixrepresentationvcsparse, +"matrixrepresentationvcsparse[op, basis] computes the matrix representation +of the operator op in a given basis (occupation number representation). +A sparse matrix is returned.", +"matrixrepresentationvc"]; +UsageWithMore[matrixrepresentationop, +"matrixrepresentationop[op, basis] computes the matrix representation +of the operator op in a given basis (creation operator representation).", +"matrixrepresentationvc"]; +UsageWithMore[makematricesbzop, +"makematricesbzop[op, bz] produces a table of operator matrices in all +invariant subspaces of the basis bz (creation operator representation)."]; +UsageWithMore[makematricesbzvc, +"makematricesbzvc[op, bz] produces a table of operator matrices in all +invariant subspaces of the basis bz (occupation number representation).", +"makematricesbzop"]; +UsageWithMore[bzvc2bzop, +"Transform a basis in occupation number representation +to a basis in creation operator representation."]; +UsageWithMore[bzop2bzvc, +"Transform a basis in creation operator representation +to a basis in occupation number representation.", "bzvc2vzop"]; +UsageWithMore[qszbasis, +"qszbasis[{ops}] returns the basis with well defined charge and +spin projection quantum numbers (Q,S_z) in creation operator representation."]; +UsageWithMore[qszbasisvc, +"qszbasisvc[{ops}] returns the basis with well defined charge and +spin projection quantum numbers (Q,S_z) in occupation number representation.", +"qszbasis"]; +UsageWithMore[qsbasis, +"qsbasis[{ops}] returns the basis with well defined charge and +total spin quantum numbers (Q,S) in creation operator representation."]; +UsageWithMore[qsbasisvc, +"qsbasisvc[{ops}] returns the basis with well defined charge and +total spin quantum numbers (Q,S) in occupation number representation."]; +UsageWithMore[zeroonvac, +"zeroonvac[expr] drops vacuum-annihilating parts of expression expr."]; +UsageWithMore[VACUUM, +"VACUUM is a placeholder for vacuum in operator expressions. If an +annihilation operator stands just before VACUUM, the expression is zero. +Likewise, if a creation operators stands just after conj[VACUUM], the +expression drops. VACUUM is normalized to 1, i.e. +vev[nc[conj[VACUUM], VACUUM]] = 1."]; +UsageWithMore[projector, +"projector[op[j], PROJ] returns the projection operator corresponding +to the fermionic operator op[j]. PROJ = PROJ0 | PROJUP | PROJDO +| PROJ2 | PROJ1 | PROJ02."]; +UsageWithMore[projector0, +"projector0[op[j]] returns the projection operator corresponding to the +fermionic operator op[j]. It projects to the subspace with zero +occupancy.", "projector"]; +UsageWithMore[projectorUP, +"projectorUP[op[j]] returns the projection operator corresponding to the +fermionic operator op[j]. It projects to the subspace with a single +spin up particle.", "projector"]; +UsageWithMore[projectorDO, +"projectorDO[op[j]] returns the projection operator corresponding to the +fermionic operator op[j]. It projects to the subspace with a single +spin down particle.", "projector"]; +UsageWithMore[projector2, +"projector2[op[j]] returns the projection operator corresponding to the +fermionic operator op[j]. It projects to the subspace with double +occupancy.", "projector"]; +UsageWithMore[projector1, +"projector1[op[j]] returns the projection operator corresponding to the +fermionic operator op[j]. It projects to the subspace with single +occupancy.", "projector"]; +UsageWithMore[projector02, +"projector02[op[j]] returns the projection operator corresponding to the +fermionic operator op[j]. It projects to the subspace with either zero +or double occupancy.", "projector"]; +UsageWithMore[decomposevc, +"decomposevc[vec, basis] decomposes a vector into components with +respect to the given basis in the occupation number representation."]; +UsageWithMore[decomposeop, +"decomposeop[vec, basis] decomposes a vector into components with +respect to the given basis in the creation operator representation.", +"decomposevc"]; +UsageWithMore[orthogvc, +"orthogvc[vecs, basis] orthogonalizes a list of vectors vecs in the +occupation number representation basis 'basis'."]; +UsageWithMore[orthogop, +"orthogop[vecs, basis] orthogonalizes a list of vectors vecs in the +creation operator representation basis 'basis'.", "orthogvc"]; +UsageWithMore[transformQStoIS, +"transformQStoIS[basis] transforms from a basis with well defined (Q,S) +quantum numbers to a basis with well defined (I,S) quantum numbers."]; +UsageWithMore[bra, +"bra[i] denotes the Dirac bra .", "bra"]; +UsageWithMore[braketrule, +"braketrule[bra[i], ket[j]] defines the value of the braket ."]; +UsageWithMore[spinbra, +"spinbra[S] returns the spin basis bras for spin S."]; +UsageWithMore[spinket, +"spinket[S] returns the spin basis kets for spin S.", "spinbra"]; +UsageWithMore[spinketbraX, +"spinketbraX[S] returns the X-component of the spin operator for spin S +in the bra-ket representation."]; +UsageWithMore[spinketbraY, +"spinketbraY[S] returns the Y-component of the spin operator for spin S +in the bra-ket representation.", "spinketbraX"]; +UsageWithMore[spinketbraZ, +"spinketbraZ[S] returns the Z-component of the spin operator for spin S +in the bra-ket representation.", "spinketbraX"]; +UsageWithMore[spinketbraP, +"spinketbraP[S] returns the spin lowering operator for spin S +in the bra-ket representation.", "spinketbraX"]; +UsageWithMore[spinketbraM, +"spinketbraM[S] returns the spin raising operator for spin S +in the bra-ket representation.", "spinketbraX"]; +UsageWithMore[spinbasis, +"spinbasis[S] returns the basis for spin S in the ket representation.", "spinbra"]; +UsageWithMore[transformfunc, +"transformfunc[list, rule] applies a rule to the states 'list' in the +occupation number representation and orthogonalizes the result."]; +UsageWithMore[transformbasis, +"transformbasis[basis, rule] applies a rule to the states in the full +basis (occupation number representation) and orthogonalizes the result."]; +UsageWithMore[mergebasis, +"mergebasis[{basis1, basis2,...}] merges several sets of basis states."]; +UsageWithMore[applybasis, +"applybasis[b, fn] applies function fn to each basis state in the +basis b."]; +UsageWithMore[transformtoPH, +"transformtoPH[b, Nph] generates a direct product of a basis b for +fermions with a basis with up-to Nph excited local phonons, i.e. Nph +is the phonon number cutoff."]; +UsageWithMore[phononnumber, +"phononnumber[Nph] returns the phonon number operator in the bra-ket +representation. Nph is the phonon number cutoff."]; +UsageWithMore[phononplus, +"phononplus[Nph] returns the phonon raising operator in the bra-ket +representation. Nph is the phonon number cutoff."]; +UsageWithMore[phononminus, +"phononminus[Nph] returns the phonon lowering operator in the bra-ket +representation. Nph is the phonon number cutoff.", "phononplus"]; +UsageWithMore[phononx, +"phononx[Nph] returns the phonon displacement operator (a^dag+a) +in the bra-ket representation. Nph is the phonon number cutoff.", +"phononplus"]; +UsageWithMore[transformtoLR, +"transformtoLR[b, map] transforms the basis b to a basis with well +defined parity. map is the list of the operators so that the sites +are symmetric with respect to the middle of this list."]; +UsageWithMore[snegold2newrules, +"Given a set of old operators o, a set of new operators n and +a set of operator expressions for the new basis, snegold2newrules[o, n, r] +produces the transformation rules from the old basis to the new basis."]; +UsageWithMore[wick, +"wick[expr] returns an operator string expressed using normal ordered +strings and contractions. This is an application of Wick's theorem."]; +UsageWithMore[contraction, +"contraction[op1, op2] contracts op1 and op2 according to the +standard definition."]; +UsageWithMore[vevwick, +"vevwick[expr] calculates the vacuum expectation value of an operator +expression using Wick's theorem.", "vev"]; +UsageWithMore[isozsq, +"isoqsz[op[j]] gives the z-component of the on-site isospin operator squared. +This term appears in the expression for e-e repulsion in the Anderson model.", +"hubbard"]; +UsageWithMore[spinspinsymmetric, +"spinspinsymmetric[{a,b,...}] returns the symmetrized spin-spin +scalar product 1/N Sum(S_a . S_b) between pairs of S=1/2 spin +operators corresponding to the levels described by the opreators +a,b,... Here N is the number of pairs."]; +UsageWithMore[maskOp, +"maskOp[op, args] calls function op while all constants declared to +sneg are masked. This prevents some interferences between sneg and +certain Mathematica functions."]; +UsageWithMore[quickISObasis, +"quickISObasis[basisops] generates an (I,S) basis."]; +UsageWithMore[majoranaQ, +"majoranaQ[z] returns True if z is a Majorana fermionic operator."]; +UsageWithMore[isnumericQ, +"isnumericQ[x] returns True if x is a numeric quantity which may +be factored out in front of a nc multiplication."]; +UsageWithMore[mcross, +"mcross[a, b] returns a cross product a x b. nc multiplication is used."]; +UsageWithMore[normvc, +"normvc[v] returns the norm of the state v in the occupation number +representation."]; +UsageWithMore[normop, +"normop[v] returns the norm of the state v in the creation operator +representation.", "normvc"]; +UsageWithMore[SimplifyKD, +"SimplifyKD[expr] simplifies KroneckerDelta and UnitStep expressions +in expr."]; +UsageWithMore[snegAssuming, +"snegAssuming[assumption, expr] evaluates expr assuming that +'assumption' holds."]; +UsageWithMore[snegSeries, +"snegSeries[f, arg, n] expands function f into a series of order n +around 0. The argument of the expansion is arg and its powers are +computed using pow."]; +UsageWithMore[SnegSimplify, +"SnegSimplify[expr] attempts to rewrite an expression in terms +of higher level sneg functions."]; +UsageWithMore[isannihilation, +"isannihilation[op] returns True if op is an annihilation operator."]; +UsageWithMore[iscreation, +"iscreation[op] returns True if op is a creation operator.", +"isannihilation"]; +UsageWithMore[dropemptysubspaces, +"dropemptysubspaces[bz] removes subspaces that contain no states +from the basis bz."]; +UsageWithMore[lrmap, +"lrmap[map] returns the substitution rules with perform a left-right +reflection transformation.", +"transformtoLR"]; +UsageWithMore[makeallmatricesbzvc, +"makeallmatricesbzvc[H, basis] computes the matrix representation +of operator H in the given basis (in the occupation number representation). +All pairs of invariant subspaces are considered."]; +UsageWithMore[makeallmatricesbzop, +"makeallmatricesbzop[H, basis] computes the matrix representation +of operator H in the given basis (in the creation operator representation). +All pairs of invariant subspaces are considered.", +"makeallmatricesbzvc"]; +UsageWithMore[asciitosneg, +"asciitosneg[string] parses a string containing a condensed operator-expression +representation and converts it in a SNEG operator-expression representation. +Operators are entered as charactes, followed by a + sign for creation +operators, the indexes are surrounded by parenthesis. The non-commutative +multiplication is implied. Bras and kets are also supported as strings. The inverse transformation can be perofmed using snegtoascii[]."]; +UsageWithMore[snegtoascii, +"snegtoascii[expr] converts a SNEG expression with non-commutative multiplication +nc[], operators, bras and kets into an ASCII string representation. If the +full conversion is not possible, the function will return a SNEG expression +featuring String objects for those parts of the expression where the conversion +was possible. The inverse transformation can be performed using asciitosneg[]."]; + +Begin["`Private`"]; (* Begin Private Context *) + +(* Sneg version *) +$SnegVersion = 1.250; + + +(** Load functions from other packages **) +(* We need KSubsets from Combinatorica by Sriram V. Pemmaraju +and Steven S. Skiena. *) +KS = Compile[ +{ + {n, _Integer}, + {k, _Integer} +}, + Module[ + { + h, + ss = Range[k], + x + }, + Table[( + h = Length[ss]; x = n; + While[x === ss[[h]], h--; x--]; + ss = Join[Take[ss, h - 1],Range[ss[[h]]+1, ss[[h]]+Length[ss]-h+1]] + ), + {Binomial[n, k]-1} + ](* end Table *) + ](* end Module *) +](* end Compile *); + +KSubsets[l_List,0] := { {} } +KSubsets[l_List,1] := Partition[l,1] +KSubsets[l_List,2] := Flatten[Table[{l[[i]], l[[j]]}, + {i, Length[l]-1}, + {j, i+1, Length[l]} + ], + 1 + ] +KSubsets[l_List,k_Integer?Positive] := {l} /; (k == Length[l]) +KSubsets[l_List,k_Integer?Positive] := {} /; (k > Length[l]) +KSubsets[s_List, k_Integer] := Prepend[Map[s[[#]] &, KS[Length[s], k]], s[[Range[k] ]] ] + +(* Mathematica >6.x required *) +(* OLD 1: snegorthog[m_] := Orthogonalize[m, Method->"Householder"]; *) +(* OLD 2: snegorthog[m_] := Orthogonalize[m, Dot[Conjugate[#1], #2] &]; *) +snegorthog[m_] := Orthogonalize[m, Simplify[Dot[Conjugate[#1], #2]] &]; +snegorthog[{{0}}] := {}; +snegorthog[m_] /; Tr[Abs[m]] == 0 := {}; + +(* When PrettyOutput is set to True, nc[] is formated in infix notation with +CenterDot and operators are formated in familiar notation with daggers, +arrows for spin projection, etc. *) +(* When PrettyInput is set to True, some input expressions are parsed, +providing a counterpart to PrettyOutput. *) + +If[!ValueQ[PrettyOutput], PrettyOutput = True]; +If[!ValueQ[PrettyInput], PrettyInput = True]; + + + +(* Standard shorthands for operator indexes; some functions depend +explicitly upon these particular values. This holds notably for ap[]. +For this reasion, these values have Protected attribute. *) + +ClearAttributes[{CR, AN, DO, UP}, Protected]; +CR = 0; +AN = 1; +DO = 0; +UP = 1; +SetAttributes[{CR, AN, DO, UP}, Protected]; + +(* list?? store symbols that have been promoted to operators, +constants or indexes. *) +If[!ValueQ[listfermionoperators], listfermionoperators={} ]; +If[!ValueQ[listbosonoperators], listbosonoperators={} ]; +If[!ValueQ[listmajoranaoperators], listmajoranaoperators={} ]; +If[!ValueQ[listspinoperators], listspinoperators={} ]; + +If[!ValueQ[listrealconstants], listrealconstants={} ]; +If[!ValueQ[listintegerconstants], listintegerconstants={} ]; +If[!ValueQ[listcomplexconstants], listcomplexconstants={} ]; +If[!ValueQ[listgrassmanconstants], listgrassmanconstants={} ]; + +If[!ValueQ[listfreeindexes], listfreeindexes={} ]; + +(* addto[list, a] adds a to list and removes duplicates *) +addto[list_, element_] := (list = Union[Append[list, element]]); +SetAttributes[addto, HoldFirst]; + +(* Define formatting rules for an operator (pretty printing). *) +SnegPPoperator[x_, spin_:1/2] := Block[{}, + (* With spin index *) + Format[x[i_, j___, sigma_], StandardForm] := + Module[{dag, sgm}, + dag = Which[i === CR, "\[Dagger]", + i === AN, "", + True, Null]; + sgm = Which[ + spin == 1/2 && sigma === UP, + (* XXX: TEMPORARY FIX: Should be \[UpArrow] *) + StyleForm["\[DoubleUpArrow]", FontColor -> Red], + spin == 1/2 && sigma === DO, + StyleForm["\[DownArrow]", FontColor -> Blue], + spin == 0 && sigma === 0, + "", + True, + sigma ]; + + DisplayForm @ If[dag === Null, + RowBox[{x, "[", RowBox[{i,j,sigma}], "]"}], + StyleBox[ Subsuperscript[x, RowBox @ Append[{j}, sgm], dag ], + ScriptSizeMultipliers -> 1, + ScriptBaselineShifts -> {1, 1}] + ] + ]; + + (* Without any indexes (apart from creation/annihilation type index) *) + Format[x[i_], StandardForm] := + Module[{dag}, + dag = Which[i === CR, "\[Dagger]", + i === AN, "", + True, Null]; + + DisplayForm @ If[dag === Null, + RowBox[{x, "[", RowBox[{i}], "]"}], + StyleBox[ Superscript[x, dag ], + ScriptSizeMultipliers -> 1, + ScriptBaselineShifts -> {1, 1}] + ] + ]; + + (* Avoid superfluous parenthesis *) + (* MakeBoxes[a_ * x[i_, j___, sigma_], fmt:(StandardForm|TraditionalForm)] := + RowBox[{MakeBoxes[a, fmt], MakeBoxes[x[i, j, sigma], fmt]}]; *) + ]; + +(* Convert x to a string in the TeX format *) +TeXString[x_] := ToString[x, TeXForm]; +SetAttributes[TeXString, Listable]; + +(* TeX conversion rules for a (fermionic) operator *) +(* XXX: broken as of 24.8.2010 due to a bug in Mathematica which reapplies TeXForm[] + to the output from our definition of TeXForm formatting. Weird! *) +SnegTeXForm[x_, spin_] := Module[{}, + Format[x[i_, j___, sigma_], TeXForm] := SnegTeXForm1[x[i, j, sigma]]; + + SnegTeXForm1[op_[i_, j___, sigma_]] := + Module[{dag, sgm}, + dag = Which[i === CR, "\\dag", + i === AN, "{}", + True, "{}"]; + sgm = Which[ + spin == 1/2 && sigma === UP, + "\\uparrow", + spin == 1/2 && sigma === DO, + "\\downarrow", + True, + sigma ]; + SequenceForm[ ToString[op], "^", dag, "_{", Append[TeXString[{j}], sgm], "}" ] + ]; + + Format[x[i_], TeXForm] := SnegTeXForm2[x[i]]; + + SnegTeXForm2[op_[i_]] := + Module[{dag}, + dag = Which[i === CR, "\\dag", + i === AN, "{}", + True, "{}"]; + SequenceForm[ ToString[op], HoldForm["^"], dag ] + ]; +]; + +(* Default spin is 1/2 *) +(* Note that the spin of a given operator needs to be redefined + by changing the appropriate upvalue, e.g. spinup[c] ^= 3/2. *) +spinof[x_] := 1/2; + +(* Add one fermion operator with given spin. *) +snegaddop[{x_, spin_}] := Block[{}, + operatorQ[x] ^= True; + fermionQ[x] ^= True; + spinof[x] ^= spin; + SetAttributes[x, NHoldAll]; + +(* SnegTeXForm[x, spin]; *) + If[PrettyOutput, SnegPPoperator[x, spin] ]; + + addto[listfermionoperators, x]; +]; + +(* By default, all fermionic operators are assumed to +be spin-1/2 operators. *) +snegaddop[x_] := snegaddop[{x, 1/2}]; + +(* Define which symbols denote fermionic operators *) +snegfermionoperators[l__] := Scan[snegaddop, {l}]; + +(* Spinless operators don't carry a spin index. This affects + their pretty printing. *) +snegspinlessfermionoperators[l__] := Scan[snegaddop[{#, 0}]&, {l}]; + +(* Define formatting rules for a bosonic operator. *) +SnegPPbosonoperator[x_] := Block[{}, + Format[x[i_, j__], StandardForm] := + Module[{dag}, + dag = Which[i === CR, "\[Dagger]", + i === AN, "", + True, Null]; + + DisplayForm @ If[dag === Null, + RowBox[{x, "[", RowBox[{i,j}], "]"}], + StyleBox[ Subsuperscript[x, RowBox @ {j}, dag ], + ScriptSizeMultipliers -> 1, + ScriptBaselineShifts -> {1, 1}] + ] + ]; + + Format[x[i_], StandardForm] := + Module[{dag}, + dag = Which[i === CR, "\[Dagger]", + i === AN, "", + True, Null]; + + DisplayForm @ If[dag === Null, + RowBox[{x, "[", RowBox[{i}], "]"}], + StyleBox[ Superscript[x, dag], + ScriptSizeMultipliers -> 1, + ScriptBaselineShifts -> {1, 1}] + ] + ]; + + (* Avoid superfluous parenthesis *) + (* MakeBoxes[a_ * x[i_, j___], fmt:(StandardForm|TraditionalForm)] := + RowBox[{MakeBoxes[a, fmt], MakeBoxes[x[i, j, sigma], fmt]}]; *) + ]; + + +(* Add one bosonic operator *) +snegaddbosonop[x_] := Block[{}, + operatorQ[x] ^= True; + bosonQ[x] ^= True; + SetAttributes[x, NHoldAll]; + + If[PrettyOutput, SnegPPbosonoperator[x] ]; + + addto[listbosonoperators, x]; +]; + +(* Define which symbols denote bosonic operators *) +snegbosonoperators[l__] := Scan[snegaddbosonop, {l}]; + +(* Define which symbols are numerical constants *) +snegrealconstants[l__] := Scan[ +{ (* We must define the symbol to be numeric. *) + (* Setting NumberQ[#] ^= True; is not sufficient! *) + NumericQ[#] ^= True; + isnumericQ[#] ^= True; + Conjugate[#] ^= #; + addto[listrealconstants, #]; + + Unprotect[Power]; + Power[Power[#, n_], 1/2] /; (EvenQ[n]) := Power[#,n/2]; + Power[-Power[#, n_], 1/2] /; (EvenQ[n]) := I Power[#,n/2]; + Protect[Power]; +}&, {l}]; + +sneggrassmanconstants[l__] := Scan[ +{ + grassmanQ[#] ^= True; + grassmanQ[conj[#]] ^= True; + If[PrettyOutput, + Format[conj[#]] := OverBar[#]; + ]; + addto[listgrassmanconstants, #]; +}&, {l}]; + +snegmajoranaoperators[l__] := Scan[ +{ + operatorQ[#] ^= True; + majoranaQ[#] ^= True; + conj[#] ^= #; + SetAttributes[#, NHoldAll]; + addto[listmajoranaoperators, #]; +}&, {l}]; + + +(* CONVENTION: the last index is {x=1,y=2,z=3,p=4,m=5}. *) +SPININDEXx=1; +SPININDEXy=2; +SPININDEXz=3; +SPININDEXp=4; +SPININDEXm=5; + +(* To save some typing.. One must be careful not to redefine x,y,z,p and m. *) +(* XXX: Protect[] these symbols?? *) +(* Commented out, Nov 2012 *) +(* +op_?spinQ[i___, x] ^= op[i, SPININDEXx]; +op_?spinQ[i___, y] ^= op[i, SPININDEXy]; +op_?spinQ[i___, z] ^= op[i, SPININDEXz]; +op_?spinQ[i___, p] ^= op[i, SPININDEXp]; +op_?spinQ[i___, m] ^= op[i, SPININDEXm]; +*) + +(* Define formatting rules for a spin operator. *) +SnegPPspinoperator[x_] := Block[{}, + Format[x[i___, j_], StandardForm] := + Module[{ndx}, + ndx = Which[j === SPININDEXx, StyleForm["x", FontColor -> Blue], + j === SPININDEXy, StyleForm["y", FontColor -> Blue], + j === SPININDEXz, StyleForm["z", FontColor -> Blue], + j === SPININDEXp, StyleForm["+", FontColor -> Blue], + j === SPININDEXm, StyleForm["-", FontColor -> Blue], + True, j]; + + DisplayForm @ StyleBox[ + If[ {i} === {}, Superscript[x, ndx], + Subsuperscript[x, RowBox @ {i}, ndx ] ], + ScriptSizeMultipliers -> 1, + ScriptBaselineShifts -> {1, 1}] + ]; + + + (* Avoid superfluous parenthesis *) + (* MakeBoxes[a_ * x[i_, j___], fmt:(StandardForm|TraditionalForm)] := + RowBox[{MakeBoxes[a, fmt], MakeBoxes[x[i, j], fmt]}]; *) + ]; + +snegspinoperators[l__] := Scan[ +{ + operatorQ[#] ^= True; + spinQ[#] ^= True; + SetAttributes[#, NHoldAll]; + + (* This allows for some tricks later on... *) + #[i___, -SPININDEXx] = -#[i, SPININDEXx]; + #[i___, -SPININDEXy] = -#[i, SPININDEXy]; + #[i___, -SPININDEXz] = -#[i, SPININDEXz]; + #[i___, 0] = 0; + + If[PrettyOutput, SnegPPspinoperator[#] ]; + + addto[listspinoperators, #]; +}&, {l}]; + + +snegintegerconstants[l__] := Scan[ +{ IntegerQ[#] ^= True; + isnumericQ[#] ^= True; + Conjugate[#] ^= #; + addto[listintegerconstants, #];}&, +{l}]; + +snegcomplexconstants[l__] := Scan[ +{ NumericQ[#] ^= True; + isnumericQ[#] ^= True; + addto[listcomplexconstants, #];}&, +{l}]; + +(* Define which symbols are free indexes that appear in sum[]s *) +snegfreeindexes[l__] := Scan[ +{ NumericQ[#] ^= True; + isnumericQ[#] ^= True; + freeindex[#] ^= True; + addto[listfreeindexes, #];}&, +{l}]; + +(* Make commonly used functions Numeric *) +SetAttributes[KroneckerDelta, NumericFunction]; +SetAttributes[UnitStep, NumericFunction]; + +Unprotect[KroneckerDelta]; +KroneckerDelta[a_List, b_List] /; Length[a] == Length[b] := + Inner[KroneckerDelta, a, b, Times]; + +KroneckerDelta /: Conjugate[KroneckerDelta[a__]] := + KroneckerDelta[a] /; Conjugate[{a}] === {a}; + +Protect[KroneckerDelta]; + +Unprotect[UnitStep]; +UnitStep /: Conjugate[UnitStep[a_]] := + UnitStep[a] /; Conjugate[a] === a; +Protect[UnitStep]; + +If[PrettyOutput, + Unprotect[KroneckerDelta]; + Format[KroneckerDelta[x_, y_], StandardForm] := + DisplayForm @ + StyleBox[SubscriptBox[StyleForm["\[Delta]", FontColor -> Green], + RowBox[{x,y}]], + ScriptSizeMultipliers -> 1, + ScriptBaselineShifts -> {1, 1}]; + Protect[KroneckerDelta]; + + Unprotect[UnitStep]; + Format[UnitStep[x_], StandardForm] := DisplayForm @ + RowBox[{StyleForm["\[Theta]", FontColor -> Green], + "(", x, ")"}]; + Protect[UnitStep]; +]; + +(* Checking for operator quantities with indexes *) +operatorQ[op_?operatorQ[k___]] := True; +fermionQ[op_?fermionQ [k___]] := True; +bosonQ[op_?bosonQ [k___]] := True; +majoranaQ[op_?majoranaQ[k___]] := True; +spinQ[op_?spinQ[k___]] := True; + +(* Fall back: a fermionic operator is clearly an operator. *) +operatorQ[op_?fermionQ] := True; +operatorQ[op_?bosonQ] := True; +operatorQ[op_?majoranaQ] := True; +operatorQ[op_?spinQ] := True; + +(* isnumericQ[] return True if the expression does not involve operators. + The default behavior is to check numericity using NumericQ and NumberQ, + but it is possible to override this function for particular purposes. *) + +isnumericQ[z_] := NumericQ[z] || NumberQ[z]; + +(* Make "op" be a linear operator in all of its arguments. *) +(* Such operators are, for example, nc[], dd[], ap[], invertspin[]... *) + +sneglinearoperator[op_] := { + op[a___, z_?isnumericQ, b___] := z op[a, b]; + op[a___, z_?isnumericQ b_, c___] := z op[a, b, c]; + op[a___, b1_ + b2_, c___] := op[a, b1, c] + op[a, b2, c]; +(* The alternative version using Distribute shows no speed improvement... *) +(* op[a___, b_Plus, c___] := Distribute[Unevaluated[op[a, b, c]], Plus, op]; *) +}; + +(* Similar to sneglinearoperator[op], but the operator op is + made linear in the first argument only. Such operators are, + for example, vev[], zeroonvac[].*) + +sneglinearoperatorFirst[op_] := { + op[z_?isnumericQ] := z; + op[z_?isnumericQ, b__] := z op[b]; + op[z_?isnumericQ b_, c___] := z op[b, c]; (* Important: b_, not b__ ! *) + op[b1_ + b2_, c___] := op[b1, c] + op[b2, c]; +}; + +(* Definition for an operator with at most two arguments; the operator +is conjugate-linear with respect to the first argument and linear +with respect to the second argument. This is the usual convention +in the Dirac's bra-ket notation, but differs from the usual convention +in the mathematical literature. *) + +snegsesquilinearoperator[op_] := { + op[z_?isnumericQ a_, b_] := conj[z] op[a, b]; + op[a_, z_?isnumericQ b_] := z op[a, b]; + op[a1_ + a2_, b_] := op[a1, b] + op[a2, b]; + op[a_, b1_ + b2_] := op[a, b1] + op[a, b2]; + (* Special case *) + op[a_, 0] := 0; + op[0, b_] := 0; +}; + +(* nc is a reserved word and it is obviously never an operator. +This definition speeds up some tests. *) +operatorQ[nc] = False; + +(**** Define operations for non-commuting operators ****) + +nc[a___, HoldPattern[nc[b__]], c___] := nc[a, b, c]; (* Associativity *) + +(* An alternative approach to make nc[] associative is to set the +attributes Flat and OneIdentity. In this case, some tweakings of the +definitions of nc[] are required in order to avoid infinite recursion +loops. I've managed to implement such a version (4.7.2006) and it +works -- much slower! (A working version of that approach is stored as +sneg-Flat-OneIdentity.m in directory special/ for experimental +purposes and further testing. It even passes all the regression +tests!). It is slower by an order of magnitude, apparently because the +pattern matching process is more involved. *) + +sneglinearoperator[nc]; + +(* nc[] is Listable *) +SetAttributes[nc, Listable]; + +(* Product of zero operators equals one, see Times[]. *) +nc[] := 1; + +(* Rule for single terms, see Times[]. *) +(* This rule implies idempotency of multiplication. *) +(* NOTE: When this rule interferes with some other desired behavior, +one can use a workaround such as nc[c,]/.Null->1 *) +nc[c_] := c; + +(* Pauli rule *) +nc[a___, op_[k___], op_[k___], b___] /; fermionQ[op] := 0; + +ordering[op_?fermionQ] = EMPTY; (* Default ordering is "Empty band" *) +(* Override with ordering[op] = SEA or some other custom definition. *) + +(* In the case of empty band, we always reorder so that the creation +operators go to the left of the operator product. This is achieved by +standard ordering test OrderedQ, since the creation operators corresponds to +CR=0, while annihilation operators correspond to AN=1. As a side-effect, the +operators are also ordered in increasing order with respect to the following +operator indexes. *) + +(* +Canonical sorting: +================= +- creation operators to the left, annihilation operators to the right +- sort by index +- sort by spin index +=> for trivial vacuum (no particles), this reduces to simple sorting + by all indexes! +*) + +snegOrderedQ[x1:op_[___], x2:op_[___]] /; + (fermionQ[op] && ordering[op] === EMPTY) := OrderedQ[{x1, x2}]; + + +(* ordering=NONE rules, added 4. 3. 2010. *) +(* This disables reordering of the operators which is useful when + working with very long operator expressions. Note that certain + assumptions about the vacuum state and the normal ordering + are implied. *) + +snegOrderedQ[x1 : op_[___], x2 : op_[___]] /; + (fermionQ[op] && ordering[op] === NONE) := True; + +isannihilation[op_[t_, j___]] /; (ordering[op] == NONE) := (t == AN); +iscreation[op_[t_, j___]] /; (ordering[op] == NONE) := (t == CR); +contraction[op_[AN, j1___], op_[CR, j2___]] /; (ordering[op] == NONE) := + If[{j1} === {j2}, 1, 0, 0]; +contraction[op_[___], op_[___]] /; (ordering[op] == NONE) := 0; + +(* NEW RULES for ordering=SEA, 25. 3. 2007 *) + +snegOrderedQ[x1:op_[i1_, k1___], x2:op_[i2_, k2___]] /; + (fermionQ[op] && ordering[op] === SEA && + iscreation[x2] && isannihilation[x1]) := False; + +snegOrderedQ[x1:op_[i1_, k1___], x2:op_[i2_, k2___]] /; + (fermionQ[op] && ordering[op] === SEA && + iscreation[x1] && isannihilation[x2]) := True; + +(* If no other pattern was matched, use standard reordering rules. The +idea is that in the absence of information about the +creation/annihilation character of the operators, we may as well +reorder the expression in a form that is convenient for further +simplifications. *) + +snegOrderedQ[x1:op_[___], x2:op_[___]] /; + (fermionQ[op] && ordering[op] === SEA) := OrderedQ[{x1, x2}]; + +(* Ordering of bosonic operators *) +snegOrderedQ[x1:op_[___], x2:op_[___]] /; bosonQ[op] := + OrderedQ[{x1, x2}]; + +(* Ordering of spin operators *) +snegOrderedQ[x1:op_[___], x2:op_[___]] /; spinQ[op] := + OrderedQ[{x1, x2}]; + +(* Perform reordering for operators with equal Head. *) + +nc[a___, x1:op_[k1__], x2:op_[k2__], b___] /; + (fermionQ[op] && !snegOrderedQ[x1, x2]) := + -nc[a, x2, x1, b] + nc[a, acmt[x1, x2], b]; + +nc[a___, x1:op_[k1__], x2:op_[k2__], b___] /; + (bosonQ[op] && !snegOrderedQ[x1, x2]) := + nc[a, x2, x1, b] + nc[a, cmt[x1, x2], b]; + +(* Only reorder spin operator when the spin index is determined (x,y,z,p,m). *) +nc[a___, x1:op_[k1___, t1_], x2:op_[k2___,t2_], b___] /; + (spinQ[op] && 1 <= t1 <= 5 && 1 <= t2 <= 5 && !snegOrderedQ[x1, x2]) := + nc[a, x2, x1, b] + nc[a, cmt[x1, x2], b]; + + +(* IMPORTANT NOTE: fermionic operators with different Head are assumed +to anti-commute; bosonic operators with different Head are assumed to +commute; bosonic and fermionic operators are assumed to commute. *) + +(* Creation operators always moved to the left *) + +nc[a___, x1:op1_?fermionQ[i1_,___], x2:op2_?fermionQ[i2_,___], b___] /; + ( (op1 =!= op2) && (i1 == AN && i2 == CR) ) := -nc[a, x2, x1, b]; + +nc[a___, x1:op1_?bosonQ[i1_,___], x2:op2_?bosonQ[i2_,___], b___] /; + ( (op1 =!= op2) && (i1 == AN && i2 == CR) ) := nc[a, x2, x1, b]; + +(* In the case of the same *type* of operators, sort according to Head +(symbol), i.e. in the alphabetical order. *) + +nc[a___, x1:op1_?fermionQ[i1_,___], x2:op2_?fermionQ[i2_,___], b___] /; + ( (op1 =!= op2) && (i1 === i2) && !OrderedQ[{x1,x2}]) := + -nc[a, x2, x1, b]; + +nc[a___, x1:op1_?bosonQ[i1_,___], x2:op2_?bosonQ[i2_,___], b___] /; + ( (op1 =!= op2) && (i1 === i2) && !OrderedQ[{x1,x2}]) := + nc[a, x2, x1, b]; + + +nc[block[expr_]] := expr; + +sneginfixlist[l_List, elem_] := + Most @ Flatten[Map[{#, elem}&, l], 1]; + +If[PrettyOutput, + (* Infix notation with green center dots. *) + MakeBoxes[HoldPattern[nc[a__]], fmt:(StandardForm|TraditionalForm)] := + RowBox[sneginfixlist[ + Map[MakeBoxes[#, fmt]&, {a}], + StyleBox["\[CenterDot]", FontColor -> Green] ]]; +]; + +If[PrettyInput, + (* Interpret center dots as nc[] multiplication. *) + (* NOTE: The center dot should not be mistaken for the scalar product. *) + CenterDot[a___] := nc[a]; +]; + + +(* Bosonic operators are commuted to the left of fermionic operators. *) +nc[a___, op1_, op2_, b___] /; (fermionQ[op1] && bosonQ[op2]) := nc[a, op2, op1, b]; + + +(* Support for anti-commuting Grassman numbers *) +nc[a___, z1_?grassmanQ, z2_?grassmanQ, b___] /; (!OrderedQ[{z1,z2}]) := + -nc[a, z2, z1, b]; + +nc[a___, z_, b___, z_, c___] /; grassmanQ[z] := 0; + +nc[a___, b_, z_, c___] /; (fermionQ[b] && grassmanQ[z]) := + -nc[a, z, b, c]; + + + +(* Support for Majorana fermions *) + +nc[a___, op_, op_, c___] /; majoranaQ[op] := 1/2 nc[a, c]; + +(* bug fixed, 24.9.2013 *) +nc[a___, x1:op1_?majoranaQ, x2:op2_?majoranaQ, c___] /; + (op1 =!= op2 && !OrderedQ[{x1, x2}]) := -nc[a, op2, op1, c] + nc[a, acmt[x1, x2], c]; + +snegOrderedQ[x1:op_[__], x2:op_[__]] /; majoranaQ[op] := OrderedQ[{x1, x2}]; + +acmt[x1:op_[j1_], x2:op_[j2_]] /; majoranaQ[op] := + If[Length[{j1}] == Length[{j2}], + Inner[KroneckerDelta, {j1}, {j2}, Times], 0]; + +acmt[op1_Symbol, op2_Symbol] /; majoranaQ[op1] && majoranaQ[op2] && op1 =!= op2 := 0; + + +(* Majorana operators are always anti-commuted to the left of the (Dirac) +fermionic operators. This allows for mixed Majorana-Dirac expressions. *) +nc[a___, x1:op1_?fermionQ, x2:op2_?majoranaQ, b___] := -nc[a, x2, x1, b]; + +(* Coherent states *) +coh[z_, c_[k___]] := nc[Exp[-nc[conj[z], z]/2], Exp[nc[c[CR, k], z]], VACUUM]; + +(* Count the number of operators with anticommuting properties *) + +acmtcount[l_List] := Count[l, _?((fermionQ[#] || grassmanQ[#]) &), 1]; + + +(* Berezin integral *) + +nc[a___, int[z_], b___] /; (grassmanQ[z] && FreeQ[{b}, z, 1]) := 0; +nc[a___, int[z_], z_, b___] /; (grassmanQ[z] && FreeQ[{b}, z, 1]) := + nc[a, b]; +nc[a___, int[z_], b__, z_, c___] /; + (grassmanQ[z] && FreeQ[{b}, z, 1] && FreeQ[{c}, z, 1]) := + nc[a, b, c] (-1)^acmtcount[{b}]; + +nc[int[z__]] := 0; +nc[a___, int[z__], b___] := nc[a, Sequence @@ Map[int, {z}], b]; + + +(* NOTE: 0.->0 rules are useful to eliminate numerically zero terms in +the expressions *) +nonuls = {0. -> 0, 0. I -> 0, Complex[0.,0.] -> 0, Complex[x_, 0.]->x, + Complex[0.,y_] -> I y}; + + + +(* Baker-Hausdorff relations for simple cases *) + +(* ++ Exp Exp ++ *) + +nc[a___, Exp[b1_], Exp[b2_], c___] /; + (Simplify[ komutator[b1, b2] ] === 0) := nc[a, Exp[b1 + b2], c]; + +nc[a___, Exp[b1_], Exp[b2_], c___] /; + (Simplify[ superkomutator[b1, b2, 2] ] === 0 && + Simplify[ superkomutator[b2, b1, 2] ] === 0) := + nc[a, Exp[b1 + b2], Exp[komutator[b1, b2]/2], c]; + +(* ++ Exp op Exp ++ *) + +nc[a___, Exp[b1_], b2__, Exp[b3_], c___] /; + (b1 == -b3 && Simplify[ komutator[b1, nc[b2]] ] === 0) := + nc[a, nc[b2], c]; + +nc[a___, Exp[b1_], b2__, Exp[b3_], c___] /; + (b1 == -b3 && Simplify[ superkomutator[b1, nc[b2], 2] ] === 0) := + nc[a, nc[b2], c] + nc[a, komutator[b1, nc[b2]], c]; + +nc[a___, Exp[b1_], b2__, Exp[b3_], c___] /; (b1 == -b3) := + Module[{kmt}, + kmt = Simplify[ + (Simplify @ komutator[b1, nc[b2]] /. nonuls) / nc[b2] ]; + Exp[kmt] nc[a, nc[b2], c] /; NumericQ[kmt] + ]; + +nc[a___, Exp[b1_], b2__, Exp[b3_], c___] /; (b1 == -b3) := + Module[{kmt}, + kmt = Simplify[ + (Simplify @ superkomutator[b1, nc[b2], 2] /. nonuls) / nc[b2] ]; + nc[a, Simplify[ Cos[I Sqrt[kmt]] nc[b2] + + Sin[I Sqrt[kmt]]/(I Sqrt[kmt]) * + komutator[b1, nc[b2]]] /. nonuls, c] /; + NumericQ[kmt] + ]; + +(* ++ op Exp ++ *) + +nc[a___, b2__, Exp[b3_], c___] /; + Simplify[ komutator[b3, nc[b2]] ] === 0 := + nc[a, Exp[b3], nc[b2], c]; + +nc[a___, b2__, Exp[b3_], c___] /; + Simplify[ superkomutator[b3, nc[b2], 2] ] === 0 := + nc[a, Exp[b3], nc[b2], c] - nc[a, Exp[b3], komutator[b3, nc[b2]], c]; + +(* Mendas-Milutinovic relations, see J. Phys. A: Math. Gen. 22, L687 +(1989). *) + +nc[a___, Exp[b1_], b2__, Exp[b3_], c___] /; + (b1 == -b3 && antikomutator[b1, nc[b2]] === 0) := + nc[a, nc[b2], Exp[-2 b1], c]; + +nc[a___, Exp[b1_], b2__, Exp[b3_], c___] /; + (b1 == -b3 && superantikomutator[b1, nc[b2], 2] === 0) := + nc[a, nc[b2], Exp[-2 b1], c] + nc[a, antikomutator[b1, nc[b2]], Exp[-2 b1], c]; + +(* Special simplification rule. It is applicable, for example, +to S=1/2 spin operators. *) + +nc[a___, Exp[b_], c___] /; (NumericQ @ Simplify @ pow[b, 2]) := + Module[{t}, + t = Simplify @ pow[b, 2]; + If[t === 0, + nc[a, 1 + b, c], + nc[a, Cosh[Sqrt[t]] + b/Sqrt[t] Sinh[Sqrt[t]], c] + ] + ]; + +(* acmt[] returns 1 if one operator is creation operator, and the other +operator is the annihilation operators, while all other indexes are equal. *) + +acmt[op_[i_, ___], op_[i_, ___]] /; fermionQ[op] := 0; + +acmt[op_[i1_], op_[i2_]] /; fermionQ[op] := If[i1 =!= i2, 1, 0]; + +acmt[op_[i1_, j1___], op_[i2_, j2___]] /; fermionQ[op] := + If[i1 =!= i2 && Length[{j1}] == Length[{j2}], + Inner[KroneckerDelta, {j1}, {j2}, Times], 0]; + + +(* cmt[op1, op2] returns the commutator [op1,op2] for bosonic operators. *) +cmt[op_[i_, ___], op_[i_, ___]] /; bosonQ[op] = 0; + +cmt[op_[CR], op_[AN]] /; bosonQ[op] = -1; +cmt[op_[AN], op_[CR]] /; bosonQ[op] = +1; + +cmt[op_[CR, j1__], op_[AN, j2__]] /; bosonQ[op] := + If[Length[{j1}] == Length[{j2}], + -Inner[KroneckerDelta, {j1}, {j2}, Times], 0]; + +cmt[op_[AN, j1__], op_[CR, j2__]] /; bosonQ[op] := + If[Length[{j1}] == Length[{j2}], + +Inner[KroneckerDelta, {j1}, {j2}, Times], 0]; + + + +(* cmt[op1, op2] returns the commutator [op1,op2] for spin operators. *) + +antisym[1,2] = 3; +antisym[2,3] = 1; +antisym[3,1] = 2; +antisym[2,1] = -3; +antisym[3,2] = -1; +antisym[1,3] = -2; +antisym[i_, i_] = 0; + +cmt[op_[i___, t1_], op_[i___, t2_]] /; spinQ[op] && 1 <= t1 <= 3 && 1 <= t2 <= 3 = + I op[i, antisym[t1, t2]]; + +cmt[op_[i___, SPININDEXz], op_[i___, SPININDEXp]] /; spinQ[op] = +op[i, SPININDEXp]; +cmt[op_[i___, SPININDEXz], op_[i___, SPININDEXm]] /; spinQ[op] = -op[i, SPININDEXm]; + +cmt[op_[i___, SPININDEXp], op_[i___, SPININDEXz]] /; spinQ[op] = -op[i, SPININDEXp]; +cmt[op_[i___, SPININDEXm], op_[i___, SPININDEXz]] /; spinQ[op] = +op[i, SPININDEXm]; + +cmt[op_[i___, SPININDEXp], op_[i___, SPININDEXp]] /; spinQ[op] = 0; +cmt[op_[i___, SPININDEXm], op_[i___, SPININDEXm]] /; spinQ[op] = 0; + +cmt[op_[i___, SPININDEXp], op_[i___, SPININDEXm]] /; spinQ[op] = +2 op[i, SPININDEXz]; +cmt[op_[i___, SPININDEXm], op_[i___, SPININDEXp]] /; spinQ[op] = -2 op[i, SPININDEXz]; + +cmt[op_[i___, SPININDEXp], op_[i___, SPININDEXx]] /; spinQ[op] = 1 op[i, SPININDEXz]; +cmt[op_[i___, SPININDEXp], op_[i___, SPININDEXy]] /; spinQ[op] = I op[i, SPININDEXz]; +cmt[op_[i___, SPININDEXm], op_[i___, SPININDEXx]] /; spinQ[op] = -1 op[i, SPININDEXz]; +cmt[op_[i___, SPININDEXm], op_[i___, SPININDEXy]] /; spinQ[op] = I op[i, SPININDEXz]; + +cmt[op_[i___, SPININDEXx], op_[i___, SPININDEXp]] /; spinQ[op] = -1 op[i, SPININDEXz]; +cmt[op_[i___, SPININDEXy], op_[i___, SPININDEXp]] /; spinQ[op] = -I op[i, SPININDEXz]; +cmt[op_[i___, SPININDEXx], op_[i___, SPININDEXm]] /; spinQ[op] = 1 op[i, SPININDEXz]; +cmt[op_[i___, SPININDEXy], op_[i___, SPININDEXm]] /; spinQ[op] = -I op[i, SPININDEXz]; + +cmt[op_[i1__, j1_], op_[i2__, j2_]] := + If[Length[{i1}] == Length[{i2}], + +Inner[KroneckerDelta, {i1}, {i2}, Times], 0] cmt[op[i1, j1], op[i1, j2]]; + +cmt[op1_?spinQ[i1___], op2_?spinQ[i2___]] /; (op1 =!= op2) = 0; + +(* (Anti-)commutators *) +SetAttributes[komutator, Listable]; +komutator[a_, b_] := nc[a, b] - nc[b, a]; + +SetAttributes[antikomutator, Listable]; +antikomutator[a_, b_] := nc[a, b] + nc[b, a]; + +(* Aliases with correct spelling *) +commutator[a___] := komutator[a]; +anticommutator[a___] := antikomutator[a]; + +(* Nested (anti-)commutators *) +superkomutator[a_, b_, 0] := b; +superkomutator[a_, b_, 1] := Simplify @ komutator[a, b]; +superkomutator[a_, b_, i_Integer] /; i > 1 := + Simplify @ komutator[a, Simplify @ superkomutator[a, b, i - 1] ]; + +superantikomutator[a_, b_, 0] := b; +superantikomutator[a_, b_, 1] := Simplify @ antikomutator[a, b]; +superantikomutator[a_, b_, i_Integer] /; i > 1 := + Simplify @ antikomutator[a, Simplify @ superantikomutator[a, b, i - 1]]; + +(* Aliases with correct spelling *) +supercommutator[a___] := superkomutator[a]; +superanticommutator[a___] := superantikomutator[a]; + +(* (Manual) normal ordering *) +SetAttributes[normalorder, Listable]; +normalorder[a_] := a - vev[a]; + +SetAttributes[normalorderwick, Listable]; +normalorderwick[a_] := a - vevwick[a]; + +(* Fix parenthesis output in pretty printing. Used for +dd[], vc[], bra[], ket[]. *) +fixparens[head_] := Module[{}, + MakeBoxes[-head[i___], fmt:(StandardForm|TraditionalForm)] := + InterpretationBox[#1, #2] & @@ { + RowBox[{MakeBoxes["-", fmt], MakeBoxes[head[i], fmt]}], + -head[i] }; + MakeBoxes[a_Complex * head[i___], fmt:(StandardForm|TraditionalForm)] := + InterpretationBox[#1, #2] & @@ { + RowBox[{"(", MakeBoxes[a, fmt], ")", MakeBoxes[head[i], fmt]}], + a head[i] }; + MakeBoxes[a_?AtomQ * head[i___], fmt:(StandardForm|TraditionalForm)] := + InterpretationBox[#1, #2] & @@ { + RowBox[{MakeBoxes[a, fmt], MakeBoxes[head[i], fmt]}], + a head[i] }; +]; + +(* dd[] denotes a double dots surrounding an expression, i.e. a normal +ordered operator string. This is merely a placeholder. To explicitly +evaluate the normal ordered expression, use normalorder[] function. *) + +If[PrettyOutput, + (* Enclose the expression between two dark green double dots. *) + Format[dd[i__], StandardForm] := DisplayForm @ + RowBox[{StyleForm["\[Colon]", FontColor -> Green], + i, + StyleForm["\[Colon]", FontColor -> Green]}]; + + fixparens[dd]; +]; + +sneglinearoperator[dd]; + +(* Associativity *) +dd[a___, HoldPattern[dd[b__]], c___] := dd[a, b, c]; + +dd[] := 1; + +(* Conjugation *) +SetAttributes[conj, Listable]; +conj[z_?isnumericQ] := Conjugate[z]; +conj[z_?isnumericQ a_] := Conjugate[z] conj[a]; +conj[a_ + b_] := conj[a] + conj[b]; +conj[nc[a_, b__]] := conj[nc[b]] ~ nc ~ conj[a]; + +(* Complex (Dirac) fermion *) +conj[op_?fermionQ[i_, j___]] := op[If[i == CR, AN, CR], j]; + +(* Boson *) +conj[op_?bosonQ[i_, j___]] := op[If[i == CR, AN, CR], j]; + +(* Real (Majorana) fermion *) +conj[x:op_?majoranaQ[___]] := x; + +(* Conjugation of sums *) +conj[sum[z_, it_List]] := sum[conj[z], it]; + +(* Conjugation of exponential terms *) +conj[Exp[expr_]]:=Exp[conj[expr]]; + +(* conj threads over Rules *) +conj[Rule[a_, b_]] := Rule[conj[a], conj[b]]; + +(* NOTE: Repeated conjugation is automatically eliminated! *) (* A +continuous linear operator defined on all of a Hilbert space has a unique +adjoint operator, so this rule is usually valid. However this does not hold +for operators that are defined only on a subset of the Hilbert space and for +unbounded operators. In the case of finite-dimensional Hilbert spaces that +are typically used in sneg, there is no difference between adjoint operators +and Hermitian conjugate operators, thus there are no reasons for concern and +the repeated conjugation may be automatically eliminated. *) + +(* See also Francois Gieres, "Mathematical surprises and Dirac's formalism +in quantum mechanics", quant-ph/990069. *) + +conj[conj[a_]] := a; + +(**** Symbolic sums ****) + +(* No sum at all *) +sum[z_, {}] := z; + +(* The indexes must be sorted for the simplifications to work correctly *) +sum[z_, it_List] := sum[z, Sort[it]] /; Not[OrderedQ[it]]; + +(* Returns True if expr contains none of symbols in list it *) +FreeOfIndexQ[expr_, it_List] := And @@ Map[FreeQ[expr, #]&, it]; + +(* Return True if expr can be factored out in front of the sum[] *) +sumFactorizableQ[expr_, it_List] := FreeOfIndexQ[expr, it]; + +(* It's always safe to extract a number *) +sum[z_ a_, it_List] := z sum[a, it] /; NumberQ[z]; + +sum[z_ a_, it_List] := z sum[a, it] /; sumFactorizableQ[z, it]; + +sum[z_, it_List] := z sum[1, it] /; (sumFactorizableQ[z, it] && z != 1); + + +(* Join two Lists. Send a message if the indexes overlap. *) +snegsumJoin::overlap = "Overlapping indexes detected when joining +`` and ``."; +snegsumJoin[l1_List, l2_List] := Module[{l}, + l = Sort @ Join[l1, l2]; + If[l =!= Union[l], + Message[snegsumJoin::overlap, l1, l2]; + ]; + l +]; + +(* Automatically rename indexes in sums when name collisions occur. +The name is formed by taking the symbol name that occurs twice and +appending a number. The number is chosen such that the new name does +not yet appear in listfreeindexes. If such behavior is undesirable, it +can be disabled by setting sumAutoRename to False. *) + +sumAutoRename = True; + +ssJWRnewname[n_, allnames_List, expr_:1] := Module[{ns, i, nn}, + ns = ToString[n]; + If[$VersionNumber >= 5.1, + (* For the following line, Mathematica 5.1 is required *) + ns = StringReplace[ns, + ShortestMatch[x : (_ ~~ ___)] ~~ DigitCharacter .. ~~ EndOfString :> x]; + ]; + For[i = 1, True, i++, + nn = ToExpression[ ns <> ToString[i] ]; + If[!MemberQ[allnames, nn] && FreeQ[expr, nn], Break[]]; + ]; + (* Automatically declare the new index quantity *) + snegfreeindexes[nn]; + nn +]; + +snegsumJoinWithRenaming[a1_, a2_, it1_List, it2_List, reverse_:False] := +Module[{l, rule}, + l = it1; + rule = {}; + Scan[ If[FreeQ[l, #] && FreeQ[a1, #], + AppendTo[l, #], + (* else *) + nn = ssJWRnewname[#, Union[it1, it2], a1]; + AppendTo[l, nn]; + AppendTo[rule, # -> nn] ]&, it2]; + If[reverse == False, + sum[nc[a1, a2 /. rule], Sort @ l], + sum[nc[a2 /. rule, a1], Sort @ l] + ] +]; + +(* Product of sums *) +nc[x1___, sum[a1_, it1_List], sum[a2__, it2_List], x2___] /; (!sumAutoRename) := + nc[x1, sum[nc[a1, a2], snegsumJoin[it1,it2]], x2]; + +nc[x1___, sum[a1_, it1_List], sum[a2_, it2_List], x2___] /; sumAutoRename := + nc[x1, snegsumJoinWithRenaming[a1, a2, it1, it2], x2]; + +(* Special case for the commutators [a,b]=ab-ba where a and b are both sums. +Renaming is performed for the second argument (b) in both cases. Typically, +the resulting expressions can be simplified to a greater extent in this +approach. *) + +komutator[x1:sum[a1_, it1_List], x2:sum[a2_, it2_List]] /; sumAutoRename := + snegsumJoinWithRenaming[a1, a2, it1, it2, False] - + snegsumJoinWithRenaming[a1, a2, it1, it2, True]; + +(* Multiple sums *) +sum[sum[a_, it1_List], it2_List] := + sum[a, snegsumJoin[it1, it2]]; + +(* Thread the first argument *) +rulesumThread = { +sum[z_List, it_List] :> Map[sum[#, it]&, z] +} + +sumThread[expr_] := expr //. rulesumThread; + +(* NOTE: Mathematica Thread[] function by default threads over all +arguments! Behavior similar to that of sumThread[] would be achieved +by using Thread[expr, List, 1]. *) + +(* Simplification rules and functions for sum expressions. *) + +rulesumExpand = { +sum[a_ + b_, it_List] :> sum[a, it] + sum[b, it], + +sum[a_, it_List] :> sum[ExpandAll[a], it] +}; + +sumExpand[expr_] := expr //. rulesumExpand; + +rulesumCollect = { +Plus[z1_. sum[a1_, it_List], z2_. sum[a2_, it_List]] :> sum[z1 a1+z2 a2, it], + +Plus[z1_. sum[a1_, it1_List], z2_. sum[a2_, it2_List]] /; + Intersection[it1, it2] =!= {} :> + Module[{int}, + int = Intersection[it1, it2]; + sum[sum[z1 a1, Complement[it1, int]] + +sum[z2 a2, Complement[it2, int]], int] + ], + +nc[a___, sum[b_,it_List], c___] :> sum[nc[a, b, c], it] +}; + +sumCollect[expr_] := expr //. rulesumCollect; + +rulesumSimplifyKD = { + sum[KroneckerDelta[n1_, n2_] a_., it_List] /; MemberQ[it, n1] :> + sum[a //. {n1 -> n2}, Complement[it, {n1}]], + + sum[KroneckerDelta[n1_, n2_] a_., it_List] /; MemberQ[it, n2] :> + sum[a //. {n2 -> n1}, Complement[it, {n2}]], + + HoldPattern[sum[a1_ + a2_, x:{i___, n1_, n2_, j___}]] /; + (FreeQ[a1, n2] && FreeQ[a2, n1]) :> + sum[1, {n2}] sum[a1, {i, n1, j}] + sum[1, {n1}] sum[a2, {i, n2, j}], + + HoldPattern[sum[a_ b_., {i___, n1_, j___}]] /; + (FreeQ[a, n1] && FreeOfIndexQ[b, {i, j}]) :> + sum[a, {i, j}] sum[b, {n1}], + + HoldPattern[sum[a_ (b_ + c_), x:{i___, n1_, j___}]] /; + (FreeQ[a, n1] && FreeOfIndexQ[b, {i, j}]) :> + sum[a c, x] + sum[a, {i, j}] sum[b, {n1}] +}; + +sumSimplifyKD[expr_] := expr //. rulesumSimplifyKD; + + + +sumAbstractIndex[expr_, ndxfunc_] := Module[{ r }, + SetAttributes[ndxfunc, NumericFunction]; + r = sum[a_, it_List] :> Module[{rule, len}, + len = Length[it]; + rule = Table[it[[i]] -> ndxfunc[i], {i, len}]; + sum[a //. rule, it //. rule] + ]; + expr /. r +]; + + + +sumNameIndex[expr_, ndxfunc_, li_List] := Module[ {len}, + len = Length[li]; + expr //. Table[ndxfunc[i] -> li[[i]], {i, len}] +]; + +rulesumStripSums = { + sum[expr_, {k___}] -> expr +}; + +(* sumFullSimplifyOLD tries to simplify an expression involving symbolic +summations. The heuristic approach used turns out to work very well in +practice. *) + +sumFullSimplifyOLD[expr_] := Module[{ex}, + ex = Expand[expr]; + ex = sumCollect[ex]; + ex = Simplify[ex]; + ex = sumSimplifyKD[ex]; + ex = SimplifyKD[ex]; + ex = Simplify[ex]; + + (* We must also expand within sum[] *) + ex = sumExpand[ex]; + ex = Simplify[ex]; + ex = sumSimplify[ex]; + ex = Simplify[ex]; + + (* and then collect again *) + ex = sumCollect[ex]; + ex = Simplify[ex] +]; + +(* Expand, simplify, collect *) +sumSimplifyFnc1[expr_] := Module[{ex}, + ex = Expand[expr]; + ex = sumExpand[ex]; + ex = Simplify[ex]; + ex = sumCollect[ex] +]; + +(* Expand within, simplify, collect *) +sumSimplifyFnc2[expr_] := Module[{ex}, + ex = sumExpand[expr]; + ex = Simplify[ex]; + ex = sumCollect[ex] +]; + +(* Simplify KD *) +sumSimplifyFnc3[expr_] := Module[{ex}, + ex = sumSimplifyKD[expr]; + ex = SimplifyKD[ex] +]; + +(* Collect, simplify, simplifyKD *) +sumSimplifyFnc4[expr_] := Module[{ex}, + ex = sumCollect[expr]; + ex = Simplify[ex]; + ex = sumSimplifyKD[ex] +]; + +(* Transform a list of rules to a list of functions that apply these rules. +Such functions can be used as transformation functions in Simplify. *) + +rulestofuncs[rules_List] := Map[Function[expr, expr /. #] &, rules]; + +(* It is better to apply individual simplification rules separately, since +the rules in each block can have conflicting effects as regards the +complexity of the expression that is being simplified. *) + +sumsimplfuncs = Join[ { Automatic, SimplifyKD, sumSimplifyFnc1, + sumSimplifyFnc2, sumSimplifyFnc3, sumSimplifyFnc4 }, + Sequence @@ Map[rulestofuncs, + {rulesumExpand, rulesumCollect, rulesumSimplifyKD} ] +]; + + +sumSimplify[expr_] := Simplify[expr, + TransformationFunctions -> sumsimplfuncs]; + +(* sumFullSimplify first calls sumFullSimplifyOLD, then applies additional + rules to furher simplify the resulting expression. *) + +sumFullSimplify[expr_] := FullSimplify[expr, + TransformationFunctions -> Prepend[sumsimplfuncs, sumFullSimplifyOLD] ]; + + +(* Annihilation vs. creation operators! *) + +isannihilation[op_[t_, j___]] /; (ordering[op] == EMPTY) := (t == AN); +iscreation[op_[t_, j___]] /; (ordering[op] == EMPTY) := (t == CR); + +(* Override the following two rules if necessary! *) +isannihilation[op_[t_, k_, j___]] /; (ordering[op] == SEA) := + (t == AN && k > 0) || (t == CR && k < 0); +iscreation[op_[t_, k_, j___]] /; (ordering[op] == SEA) := + (t == CR && k > 0) || (t == AN && k < 0); + +(* For bosonic operators *) +isannihilation[op_[t_, j___]] /; bosonQ[op] := (t == AN); +iscreation[op_[t_, j___]] /; bosonQ[op] := (t == CR); + + +(* ### Vacuum expectation value ### *) + +sneglinearoperatorFirst /@ {vev, vevwick}; + +(* Single operator (not enclosed in nc[]) *) +SetAttributes[{vev, vevwick}, Listable]; +vev[_?fermionQ[__]] := 0; +vevwick[_?fermionQ[__]] := 0; + +vev[_?bosonQ[__]] := 0; + +(* Odd number of operators *) +vev[x:HoldPattern[nc[_?fermionQ[__]..]]] /; OddQ[Length[x]] := 0; +vevwick[x:HoldPattern[nc[_?fermionQ[__]..]]] /; OddQ[Length[x]] := 0; + +vev[x:HoldPattern[nc[_?bosonQ[__]..]]] /; OddQ[Length[x]] := 0; + +(* The following rule is useful in some situations, for example when +normalorder[] is used to obtain a normal ordered operator string (which may +feature a vev term), to which we then apply vev anew: this combined +operation should yield 0. *) + +vev[vev[x__]] := vev[x]; +vev[vevwick[x__]] := vevwick[x]; +vevwick[vev[x__]] := vevwick[x]; +vevwick[vevwick[x__]] := vevwick[x]; + +(* vev is applied to addends in symbolic sums *) + +vev[sum[ops_, it_List]] := sum[ vev[ops], it ]; +vevwick[sum[ops_, it_List]] := sum[ vevwick[ops], it ]; + +(* By definition, the vacuum expectation value of a normal ordered +expression is zero! *) + +vev[dd[_]] := 0; +vevwick[dd[_]] := 0; + +(* Annihilates vacuum to the right *) + +vev[ nc[___, x:op_?fermionQ[__]] ] /; isannihilation[x] := 0; +vevwick[ nc[___, x:op_?fermionQ[__]] ] /; isannihilation[x] := 0; + +vev[ nc[___, x:op_?bosonQ[AN,___]] ] := 0; + + +(* Annihilates vacuum to the left *) + +vev[ nc[x:op_?fermionQ[__], ___] ] /; iscreation[x] := 0; +vevwick[ nc[x:op_?fermionQ[__], ___] ] /; iscreation[x] := 0; + +vev[ nc[x:op_?bosonQ[CR,___]] ] := 0; + + +(* Special rules for bosons *) + + (* VEV of a string of identical bosonic operators must be zero. + The equivalent rule for fermionic operators is not necessary, due + to Pauli exclusion rule with is enforced at all times in the definition + of nc[] multiplication. *) + +vev[ HoldPattern[ nc[op_?bosonQ[___]..] ] ] := 0; + + +(* Special rules for vev in the case of ordering=SEA *) + +(* (1-n_k) |0> *) +vev[nc[a___, op_[AN, k_, j___], op_[CR, k_, j___]]] /; + (fermionQ[op] && ordering[op] === SEA) := UnitStep[k] vev[nc[a]]; + +(* <0| (1-n_k) *) +vev[nc[op_[AN, k_, j___], op_[CR, k_, j___], b___]] /; + (fermionQ[op] && ordering[op] === SEA) := UnitStep[k] vev[nc[b]]; + +(* n_k |0> *) +vev[nc[a___, op_[CR, k_, j___], op_[AN, k_, j___]]] /; + (fermionQ[op] && ordering[op] === SEA) := UnitStep[-k] vev[nc[a]]; + +(* <0| n_k *) +vev[nc[op_[CR, k_, j___], op_[AN, k_, j___], b___]] /; + (fermionQ[op] && ordering[op] === SEA) := UnitStep[-k] vev[nc[b]]; + +(* VEV of two operators with two different sets of ultimate indexes is zero. +Here we need to be extremely careful about the free indexes which can take +an arbitrary value. Note the four-argument version of the If[] statement! *) + +vev[nc[op_[i1_, k1_, j1___], op_[i2_, k2_, j2___]]] := 0 /; + fermionQ[op] && + ordering[op] === SEA && If[{j1} != {j2}, True, False, False]; + +(* Special rule. The last index is spin and we assume that the vacuum +state has well defined z-component of spin. The first index is operator +type, which can also be interpreted as the z-component of the isospin, +and we assume that the vacuum has well define particle number. The +total isospin and spin of the operator string must be equal to zero. *) + +vev[x:HoldPattern[nc[_?((fermionQ[#] && ordering[#] === SEA)&) + [_?((# === CR || # === AN)&), _, ___, _?((# === UP || # === DO)&)]..]] ] := + Module[{tmp, isospin, spin}, + tmp = Map[(# /. opp_[ii_,kk_,jj___, ss_] -> {ii, ss}) &, List @@ x]; + isospin = tmp[[All,1]]; + spin = tmp[[All,2]]; + isospin = Map[If[# == CR, 1/2, -1/2]&, isospin]; + spin = Map[If[# == UP, 1/2, -1/2]&, spin]; + (* Note: change of spin depends on the isospin of the operator! *) + spin = spin (2isospin); (* elementwise multiplication *) + + isospin = Plus @@ isospin; + spin = Plus @@ spin; + + 0 /; (isospin != 0 || spin != 0) +]; + + +(* Finally, an explicit procedure for bilinear forms in the case of SEA +ordering. This is then used for more general expression, for example in +conjunction with Wick's theorem. *) + +vev[nc[op1_[i1_, k1_, j1___], op2_[i2_, k2_, j2___]]] := + KroneckerDelta[op1, op2] KroneckerDelta[{j1}, {j2}] KroneckerDelta[k1, k2] * + (If[i1 == CR && i2 == AN, UnitStep[-k1], 0] + + If[i1 == AN && i2 == CR, UnitStep[k1], 0]) /; + fermionQ[op1] && fermionQ[op2] && + ordering[op1] === SEA && ordering[op2] === SEA; + + +(*** Inner and outer products using nc ***) + +inner[a_, b_, c___] := Inner[nc, a, b, c]; +outer[a_, b_, c___] := Outer[nc, a, b, c]; + +(* Vector-matrix-vector product using inner multiplication *) +VMV[v1_, m_, v2_] := (v1 ~ inner ~ m) ~ inner ~ v2; + +(* Direct product of two lists. May be used in in basis construction. *) +direct[l1_List, l2_List] := Flatten @ outer[l1, l2]; + +(* DEPRECATED: mdot[a_, b_] := Inner[nc, a, b]; *) + + +(*** Auxiliary functions for operator construction ***) + +(* Number operator *) +SetAttributes[number, Listable]; + +number[op_?fermionQ[j___], sigma_] := op[CR, j, sigma] ~ nc ~ op[AN, j, sigma]; +number[op_?fermionQ[j___]] /; spinof[op] == 1/2 := + number[op[j], UP] + number[op[j], DO]; +number[op_?fermionQ[j___]] /; spinof[op] != 1/2 := + Sum[number[op[j], s], {s, -spinof[op], +spinof[op]}]; + +(* Number operator for abstract function argument *) +number[fn_Function, sigma_] := fn[CR, sigma] ~ nc ~ fn[AN, sigma]; +number[fn_Function] := number[fn, UP] + number[fn, DO]; + +(* Number operator for bosonic operators *) +number[op_?bosonQ[j___]] := op[CR, j] ~ nc ~ op[AN, j]; + +number[op_?AtomQ, a___] := number[op[], a]; + +(* Number-1 squared operator, i.e. z-component of isospin squared *) +SetAttributes[isozsq, Listable]; +isozsq[op_?fermionQ[j___]] := pow[number[op[j]]-1, 2]; + +(* Hubbard's local e-e repulsion operator *) +SetAttributes[hubbard, Listable]; +hubbard[op_?fermionQ[j___]] /; spinof[op] == 1/2 := + number[op[j], UP] ~ nc ~ number[op[j], DO]; +hubbard[op_?fermionQ[j___]] /; spinof[op] != 1/2 := Module[{smin, smax}, + smin = -spinof[op]; + smax = +spinof[op]; + Sum[Sum[number[op[j], s1] ~ nc ~ number[op[j], s2], + {s1, s2+1, smax}], {s2, smin, smax}] +]; + +hubbard[op_?fermionQ, a___] /; AtomQ[op] := hubbard[op[], a]; + +(* The Hubbard model for operator op[], nn sites and parameters + U, t and eps *) +hamiltonian[Hubbard, op_[j___], nn_, U_, t_, eps_] := + Sum[eps number[op[i, j]] + U hubbard[op[i, j]], {i, 1, nn}] + + t Sum[hop[op[i, j], op[i+1, j]], {i, 1, nn-1}]; + +(* Inter-site charge-charge interaction *) +chargecharge[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] := + number[op1[j1]] ~ nc ~ number[op2[j2]]; + +(* Spin operators (Pauli matrices) for S=1/2 fermions. Note that +they must be multiplied by 1/2 to obtain true spin operators. *) +PauliX = {{0, 1}, {1, 0}}; +PauliY = {{0, -I}, {I, 0}}; +PauliZ = {{1, 0}, {0, -1}}; +PauliPlus = PauliX + I PauliY; +PauliMinus = PauliX - I PauliY; + + + +spinxyzgen[g_] := Module[{cr, an}, + cr = {g[UP], g[DO]}; + an = conj /@ cr; + 1/2 VMV[cr, #, an]& /@ {PauliX, PauliY, PauliZ} +]; + +spinxyzgen[g_, spin_?halfintegerQ] := Module[{cr, an}, + cr = Table[g[s], {s, +spin, -spin, -1}]; + an = conj /@ cr; + VMV[cr, #, an]& /@ + {spinmatrixX[spin], spinmatrixY[spin], spinmatrixZ[spin]} +]; + +SetAttributes[{spinxyz, spinall, spinss, spinx, spiny, spinz, + spinplus, spinminus}, Listable]; + +(* Returns the x, y and z components of the spin operator. *) +spinxyz[op_?fermionQ[j___]] /; spinof[op] == 1/2 := + spinxyzgen[op[CR, j, #]&]; +spinxyz[op_?fermionQ[j___]] /; spinof[op] != 1/2 := + spinxyzgen[op[CR, j, #]&, spinof[op]]; + +spinxyz[op_] /; AtomQ[op] := spinxyz[op[]]; + +(* Spin-1/2 operator for abstract function argument. *) +spinxyz[fn_Function] := spinxyzgen[fn[CR, #]&]; + +(* Returns the x, y, and z components of the spin operator, as well as +the S^2 operator *) +spinall[op_?fermionQ[j___]] := Module[{sxyz, ss}, + sxyz = spinxyz[op[j]]; + ss = sxyz ~ inner ~ sxyz; + Append[sxyz, ss] +]; + + +(* Shorthand calls for generation of individual spin operators. *) +spinss[op_?fermionQ[j___]] := spinall[op[j]] [[4]]; +spinx[op_?fermionQ[j___]] := spinall[op[j]] [[1]]; +spiny[op_?fermionQ[j___]] := spinall[op[j]] [[2]]; +spinz[op_?fermionQ[j___]] := spinall[op[j]] [[3]]; + +spinplus[op_?fermionQ[j___]] := Simplify[ spinx[op[j]] + I spiny[op[j]] ]; +spinminus[op_?fermionQ[j___]] := Simplify[ spinx[op[j]] - I spiny[op[j]] ]; + +spinss[op_] /; AtomQ[op] := spinss[op[]]; +spinx[op_] /; AtomQ[op] := spinx[op[]]; +spiny[op_] /; AtomQ[op] := spiny[op[]]; +spinz[op_] /; AtomQ[op] := spinz[op[]]; +spinplus[op_] /; AtomQ[op] := spinplus[op[]]; +spinminus[op_] /; AtomQ[op] := spinminus[op[]]; + +(* Scalar product of spin operators for op1 and op2 *) + +spinspin[fn1_Function, fn2_Function] := Module[{sp1, sp2}, + sp1 = spinxyz[fn1]; + sp2 = spinxyz[fn2]; + Simplify[ sp1 ~ inner ~ sp2 ] +]; + +spinspin[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] := + Simplify[ spinxyz[op1[j1]] ~ inner ~ spinxyz[op2[j2]] ]; + +spinspin[op1_, op2_] /; AtomQ[op1] && AtomQ[op2] := spinspin[op1[], op2[]]; + +spinspin[op1_?fermionQ[j1___], fn2_Function] := + spinspin[op1[#1, j1, #2]&, fn2]; + +spinspin[fn1_Function, op2_?fermionQ[j2___]] := + spinspin[fn1, op2[#1, j2, #2]&]; + +spinspin[op1_?spinQ[j1___], op2_?spinQ[j2___]] := + Sum[nc[op1[j1, i], op2[j2, i]], {i, 3}]; + +(* Recall: S1 . S2 = S1^z S2^z + 1/2 (S1^+ S2^- + S1^- S2^+) *) +(* Note: spinspin[a,b] = spinspinz[a,b] + spinspinxy[a,b] -> no factor 1/2 + in this definition! *) + +(* Transverse part of the spin-spin interaction *) +spinspinpm[op1_?fermionQ[j1___],op2_?fermionQ[j2___]] := + spinplus[op1[j1]] ~ nc ~ spinminus[op2[j2]]; + +spinspinmp[op1_?fermionQ[j1___],op2_?fermionQ[j2___]] := + spinminus[op1[j1]] ~ nc ~ spinplus[op2[j2]]; + +spinspinx[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] := + Expand[ spinx[op1[j1]]~nc~spinx[op2[j2]] ]; + +spinspiny[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] := + Expand[ spiny[op1[j1]]~nc~spiny[op2[j2]] ]; + +spinspinxy[op1_?fermionQ[j1___],op2_?fermionQ[j2___]] := Module[{sp1, sp2}, + sp1 = spinspinpm[op1[j1], op2[j2]]; + sp2 = spinspinmp[op1[j1], op2[j2]]; + (sp1+sp2)/2 +]; + +(* Longitudinal part of the spin-spin interaction *) +spinspinz[op1_?fermionQ[j1___],op2_?fermionQ[j2___]] := Module[{}, + spinz[op1[j1]] ~ nc ~ spinz[op2[j2]] +]; + +(* Given a list of operators 'ops', manyspin[] constructs the +components of the total spin, as well as the total spin squared operator. *) +manyspin[ops_List] := Module[{sxyz, ss}, + sxyz = Expand @ (Plus @@ Map[spinxyz, ops]); + ss = Expand @ (sxyz ~ inner ~ sxyz); + Append[sxyz, ss] +]; + +(* Returns True if z is half-spin or entire spin *) +halfintegerQ[z_] := IntegerQ[2z]; + +(* Spin matrices for arbitrary spin *) +spinmatrixZ[S_?halfintegerQ] := + Table[If[sz1 == sz2, sz1, 0], {sz1, S, -S, -1}, {sz2, S, -S, -1}]; +spinmatrixP[S_?halfintegerQ] := + Table[If[sz1 == sz2 + 1, Sqrt[S(S + 1) - sz2(sz2 + 1)], 0], {sz1, + S, -S, -1}, {sz2, S, -S, -1}]; +spinmatrixM[S_?halfintegerQ] := + Table[If[sz1 == sz2 - 1, Sqrt[S(S + 1) - sz2(sz2 - 1)], 0], {sz1, + S, -S, -1}, {sz2, S, -S, -1}]; +spinmatrixX[S_?halfintegerQ] := 1/2(spinmatrixP[S] + spinmatrixM[S]); +spinmatrixY[S_?halfintegerQ] := 1/(2I)(spinmatrixP[S] - spinmatrixM[S]); + + +(* Nambu spinor *) +SetAttributes[nambu, Listable]; +nambu[op_?fermionQ[j___], n_:0] := {op[CR, j, UP], (-1)^n op[AN, j, DO]}; +nambu[op_?AtomQ, n___] := nambu[op[], n]; + +(* X, Y and Z componentes of the isospin operator. *) + +SetAttributes[{isospinxyz, isospin, isospinx, isospiny, isospinz, + isospinplus, isospinminus}, Listable]; + +isospinxyz[op_?fermionQ[j___], n_:0] := Module[{nam, nambar}, + nam = nambu[op[j], n]; + nambar = conj /@ nam; + 1/2 VMV[nam, #, nambar]& /@ {PauliX, PauliY, PauliZ} +]; + +isospin[x__] := isospinxyz[x]; + +(* Shorthand calls for generating individual isospin operators. *) +isospinx[op_?fermionQ[j___], n_:0] := isospin[op[j], n] [[1]]; +isospiny[op_?fermionQ[j___], n_:0] := isospin[op[j], n] [[2]]; +isospinz[op_?fermionQ[j___], n_:0] := isospin[op[j], n] [[3]]; + +isospinplus[a__] := Expand[isospinx[a] + I isospiny[a]]; +isospinminus[a__] := Expand[isospinx[a] - I isospiny[a]]; + +isospinx[op_, n___] /; AtomQ[op] := isospinx[op[], n]; +isospiny[op_, n___] /; AtomQ[op] := isospiny[op[], n]; +isospinz[op_, n___] /; AtomQ[op] := isospinz[op[], n]; +isospinplus[op_, n___] /; AtomQ[op] := isospinplus[op[], n]; +isospinminus[op_, n___] /; AtomQ[op] := isospinminus[op[], n]; + +(* Sum of isospin generators for operators enlisted in ops. *) +manyisospin::unequallength = "Different lengths: ops=`` ns=``."; +manyisospin[ops_List, ns1_:Null] := Module[{len, ns, ixyz, ii}, + len = Length[ops]; + If[ns1 === Null, ns = Table[0, {len}], ns = ns1]; + If[Length[ns] != len, + Message[manyisospin::unequallength, ops, ns]; + Return[]; + ]; + ixyz = Expand @ (Plus @@ Table[isospin[ ops[[i]], ns[[i]] ], {i, len}]); + ii = (ixyz ~ inner ~ ixyz) // Expand; + Append[ixyz, ii] +]; + +(* Basis states for a single site *) +SetAttributes[basis, Listable]; +basis[op_?fermionQ[j___]] := + {1, op[CR,j,UP], op[CR,j,DO], nc[op[CR,j,DO], op[CR,j,UP]]}; + +(* Integer powers of operators *) +SetAttributes[pow, Listable]; +pow[op_, 0] := 1; +pow[op_, 1] := op; +pow[op_, i_Integer] /; i>1 := nc[pow[op, i-1], op]; + +pow::negativepower = "Negative power detected: (``)^(``)"; +pow[op_, i_] /; i<0 := (Message[pow::negativepower, op, i]; Null); + +SetAttributes[fastpow, Listable]; +fastpow[op_, 0] := 1; +fastpow[op_, 1] := Simplify[op]; +fastpow[op_, i_Integer] /; i>1 := Simplify @ nc[fastpow[op, i-1], Simplify[op]]; + +fastpow::negativepower = "Negative power detected: (``)^(``)"; +fastpow[op_, i_] /; i<0 := (Message[fastpow::negativepower, op, i]; Null); + +(* Additional optimizations for low powers. They make use of Mathematica's + result-caching capabilities. *) +fastpow[op_, 4] := Simplify @ nc[fastpow[op, 2], fastpow[op, 2]]; +fastpow[op_, 6] := Simplify @ nc[fastpow[op, 4], fastpow[op, 2]]; +fastpow[op_, 8] := Simplify @ nc[fastpow[op, 4], fastpow[op, 4]]; + +(* DANGEROUS: avoid using Power directly! *) +(* The modified Power operators is can be used only when defined explicitly + by calling the function definencPower[]. *) +definencPower[] := Module[{}, + Unprotect[Power]; + Power[expr_, i_Integer] := pow[expr, i] /; + (!isnumericQ[expr] && + (Count[expr, nc[__]] != 0 || Count[expr, _?operatorQ] != 0) && + i >= 0); + Protect[Power]; +]; + +(* Hopping operator *) +SetAttributes[hop, Listable]; +hop[op1_?fermionQ[j1___], op2_?fermionQ[j2___], sigma_] := + op1[CR, j1, sigma] ~ nc ~ op2[AN, j2, sigma] + + op2[CR, j2, sigma] ~ nc ~ op1[AN, j1, sigma]; + +hop[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] /; + (spinof[op1] == spinof[op2] == 1/2) := + hop[op1[j1], op2[j2], UP] + hop[op1[j1], op2[j2], DO]; + +hop[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] /; + (spinof[op1] == spinof[op2] != 1/2) := + Sum[ hop[op1[j1], op2[j2], s], {s, -spinof[op1], spinof[op1]} ]; + +hop[fn1_Function, fn2_Function, sigma_] := + fn1[CR, sigma] ~ nc ~ fn2[AN, sigma] + + fn2[CR, sigma] ~ nc ~ fn1[AN, sigma]; + +hop[fn1_Function, fn2_Function] := + hop[fn1, fn2, UP] + hop[fn1, fn2, DO]; + +hop[fn1_Function, op2_?fermionQ[j2___]] := hop[fn1, op2[#1, j2, #2]&]; +hop[op1_?fermionQ[j1___], fn2_Function] := hop[op1[#1, j1, #2]&, fn2]; + +(* Hopping with spin-flip *) +SetAttributes[spinfliphop, Listable]; +spinfliphop[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] := Sum[ + op1[CR, j1, sigma] ~ nc ~ op2[AN, j2, 1-sigma] + + op2[CR, j2, sigma] ~ nc ~ op1[AN, j1, 1-sigma], +{sigma, DO, UP} ]; + +(* Hopping with hole operators: +c^\dag_{k,\sigma} f_\sigma} + f^\dag_\sigma c_{k\sigma} +\to +h_{k,-\sigma} f_\sigma + f^\dag_\sigma h^\dag_{k,-\sigma} +*) + +SetAttributes[holehop, Listable]; +holehop[oph_?fermionQ[jh___], opf_?fermionQ[jf___], sigma_] := +oph[AN, jh, sigma] ~ nc ~ opf[AN, jf, 1-sigma] + +opf[CR, jf, sigma] ~ nc ~ oph[CR, jh, 1-sigma]; + +holehop[oph_?fermionQ[jh___], opf_?fermionQ[jf___]] := + Sum[holehop[oph[jh], opf[jf], sigma], {sigma, DO, UP}]; + +(* 2-electron hopping operator *) +SetAttributes[twohop, Listable]; +twohop[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] := + nc[ op1[CR, j1, DO], op1[CR, j1, UP], op2[AN, j2, DO], op2[AN, j2, UP] ] + + nc[ op2[CR, j2, DO], op2[CR, j2, UP], op1[AN, j1, DO], op1[AN, j1, UP] ]; + +(* Hopping with a phase change *) +hopphi[op1_?fermionQ[j1___], op2_?fermionQ[j2___], sigma_, phi_] := + Exp[I phi] op1[CR, j1, sigma] ~ nc ~ op2[AN, j2, sigma] + + Exp[-I phi] op2[CR, j2, sigma] ~ nc ~ op1[AN, j1, sigma]; + +hopphi[op1_?fermionQ[j1___], op2_?fermionQ[j2___], phi_] /; + (spinof[op1] == spinof[op2] == 1/2) := + hopphi[op1[j1], op2[j2], UP, phi] + + hopphi[op1[j1], op2[j2], DO, phi]; + +(* Current operator *) +(* + H = hopphi[d[1], d[2], phi] + j = d/dt N1 (i.e. current from orbital 2 to orbital 1!!) + j = i[H,N1] = current[d[1], d[2], phi] +*) + +currentphi[op1_?fermionQ[j1___], op2_?fermionQ[j2___], sigma_, phi_] := + hopphi[op1[j1], op2[j2], sigma, phi-Pi/2]; + +currentphi[op1_?fermionQ[j1___], op2_?fermionQ[j2___], phi_] /; + (spinof[op1] == spinof[op2] == 1/2) := + currentphi[op1[j1], op2[j2], UP, phi] + + currentphi[op1[j1], op2[j2], DO, phi]; + +current[op1_?fermionQ[j1___], op2_?fermionQ[j2___], sigma_] := + currentphi[op1[j1], op2[j2], sigma, 0]; + +current[op1_?fermionQ[j1___], op2_?fermionQ[j2___]] := + currentphi[op1[j1], op2[j2], 0]; + +(****************** MISC ********************) + +sneglinearoperator[invertspin]; +SetAttributes[invertspin, Listable]; + +invertspin[] := 1; +invertspinindex[j_] := (1-j); (* This works for S=1/2 only! *) + +(* Inverts the spin of all operators (spin reflection) *) +invertspin[op_?fermionQ[i_, k___, j_]] := + op[i, k, invertspinindex[j]]; + +invertspin[HoldPattern[ nc[a__] ]] := nc @@ Map[invertspin, {a}]; + + + +(************** AP **************) + +(* mbfunc[] expands each generic operator a[j] into a pair of creation +operators for both spin orientations. This function is used in basis +construction in makebasis[] and in lrmap[] for constructing a mapping +between left and right sites under a mirror symmetry. It can be +generalized for other symmetries (i.e. orbital indexes) and other +representations of the SU(2) group (i.e. S=1 representation instead of +S=1/2). *) + +SetAttributes[mbfunc, Listable]; +mbfunc[op_[j___]] /; spinof[op] == 1/2 := + { op[CR, j, UP], op[CR, j, DO] }; +mbfunc[op_[j___]] /; spinof[op] != 1/2 := + Table[ op[CR, j, s], {s, +spinof[op], -spinof[op], -1} ]; + + +(* BASIS holds a (global) list of creation operators for all possible + single-particle states. It is constructed using a call to makebasis[], + see below. Creation operators for all many-particle states in the (Q,S) + basis are constructed by combining these operators using suitable + Clebsch-Gordan coefficients. *) + +If[!ValueQ[BASIS], BASIS = {}]; + +(* Construct basis for common enlisted operators. + Example: makebasis[{a[],f[1]}] + Expansion over both spin orientations is automatically done. *) + + +makebasis[l_List] := BASIS = Flatten @ (mbfunc @ l); +makebasis[l__] := makebasis[{l}]; + +(* Is the operator in the basis? Disregards CR vs. AN type. *) +(* NOTE: op_fermionQ would be inappropriate here! *) +goodopQ[op_[i_, j___]] := Count[BASIS, op[CR, j]] == 1; + +(* Mapping of operators to positions in state vectors. An index + is only returned if the operator actually exists in the basis, + otherwise the function call is unevaluated. *) +SetAttributes[{op2ndx, ndx2op}, Listable]; +op2ndx[op_?fermionQ[i_, j___]] := Position[BASIS, op[CR, j]] [[1, 1]] /; + goodopQ[op[i, j]]; + +ndx2op[n_Integer] := BASIS[[n]]; + + +(* Vector notation *) + +(* The default coloring function for vector objects assumes that the +basis is formed by spin-1/2 fermions and eventually an additional +element from the ket basis. *) + +vccolorDefault[i_List] := Module[{a, b}, + a = If[OddQ[Length[i]], Most[i], i]; + a = a /. {0 -> "\[EmptyRectangle]", 1 -> "\[FilledRectangle]"}; + b = Flatten[Map[{StyleForm[First[#], FontColor -> Red], + StyleForm[Last[#], FontColor -> Blue]} &, Partition[a, 2]], 1]; + If[OddQ[Length[i]], AppendTo[b, Last[i]]]; + b +]; + +(* Empty or filled black rectangles. *) + +vccolorSimple[i_List] := + i /. {0 -> "\[EmptyRectangle]", 1 -> "\[FilledRectangle]"}; + +vccolor := vccolorDefault; + +If[PrettyOutput, + Unprotect[vc]; + + Format[vc[i__], StandardForm] := + DisplayForm @ RowBox[Join[ + {StyleForm["\[LeftDoubleBracketingBar]", FontColor -> Green]}, + vccolor[{i}], + {StyleForm["\[RightAngleBracket]", FontColor -> Green]}] + ]; + + fixparens[vc]; + + Protect[vc]; +]; + +(* Prevent overwriting! *) +Protect[vc]; + +(* Vacuum vector: vc[0, 0, ..., 0] *) +vacuum[] := vc @@ Table[0, {Length[BASIS]}]; + +(* Special rule for computing scalar products *) +nc[a___, conj[v1_vc], v2_vc, b___] := scalarproductvc[v1, v2] nc[a, b]; + +(*** Applies an operator to a vector in occupation number representation ***) + +sneglinearoperator[ap]; +SetAttributes[ap, Listable]; + +ap[v_vc] := v; + +ap[HoldPattern[x:nc[___]], v_vc] := + Fold[ ap[#2, #1]&, v, Reverse[ List @@ x] ]; + + +(* NOTE: the following implementation of ap[] is valid for any spin of +the fermionic operators !!! *) + +ap[op_?fermionQ[i_, j___], v_vc] := Module[{ndx, sign, ok}, + ndx = op2ndx[op[i, j]]; + sign = If[EvenQ[Length[Cases[Take[v, ndx-1], 1]]], 1, -1]; + ok = If[v[[ndx]] == i, 1, 0]; + ok sign MapAt[(1-#)&, v, ndx] +] /; goodopQ[op[i, j]]; + +(* Resolve unevaluated repeated applications of ap by multiplying + the operators. This is particularly useful for applying Dirac bras. + See below. *) +ap[a___, HoldPattern[ ap[b___, c_] ]] := ap[nc[a,b], c]; + +(* Special rules *) +ap[1, VACUUM] := vacuum[]; +ap[a__, VACUUM] := ap[a, vacuum[]]; + +(* If a vc vector appears as the last element in a nc product, ap[] is + called automatically. This allows for equivalent treatment of basis + sets in occupation number and in creation operator representation + as regards the application of operators expression to the basis + states. *) +nc[a___, v_vc] := ap[nc[a], v]; + +(* vc2ops[] transforms a vc[] expression (occupation number +representation) for a state into an equivalent expression in terms of +the creation operators. *) + +sneglinearoperator[vc2ops]; +SetAttributes[vc2ops, Listable]; + +(* Special case: expressions with a trailing ket vector. *) +vc2ops[vc[l__, ket[k__]]] := nc[vc2ops[vc[l]], ket[k]]; + +vc2ops::incom = "Incompatible lenghts of vc[] and BASIS."; +vc2ops[x:vc[l__]] := If[Length[{l}] != Length[BASIS], + Message[vc2ops::incom]; Infinity, + nc @@ Pick[BASIS, x, 1] + ]; + +(* Inverse operation *) +ops2vc[op_, vak_] := ap[op, vak]; +ops2vc[op_] := ap[op, vacuum[]]; + +bzvc2bzop[bvc_] := Map[{First[#], vc2ops[#[[2]]]}&, bvc]; +(* If the vacuum state vac is not provided, vacuum[] is used. *) +bzop2bzvc[bz_, vak_] := Map[{First[#], ap[#[[2]], vak]}&, bz]; +bzop2bzvc[bz_] := bzop2bzvc[bz, vacuum[] ]; + + +(* It is very important that the integers in vc[] are not converted to real +values by application of N[] to some state. *) +SetAttributes[vc, NHoldAll]; + +(* Scalar product of an expressions with number representation vectors *) +(* IMPORTANT: by definition, an inner product is linear in the first argument +and conjugate-linear in the second argument, i.e., it is a sesquilinear form! +In physics, on the other hand, the convention is that the first argument is +conjugate-linear (bra), while the second argument is linear (ket). +As of sneg.m v202 (17.11.2009), scalarproductvc[] follows the second convention, +in order to be in conformance with the definition of scalarproductop[]. *) + +snegsesquilinearoperator[scalarproduct]; +snegsesquilinearoperator[scalarproductvc]; +SetAttributes[{scalarproduct, scalarproductvc, scalarproductop}, Listable]; + +scalarproductvc[a_vc, b_vc] := If[a === b, 1, 0]; + +scalarproductvc[a_, 0.] := 0; +scalarproductvc[0., b_] := 0; +scalarproductvc[a_, 0.+I 0.] := 0; +scalarproductvc[0.+I 0., b_] := 0; + +scalarproduct[a___] = scalarproductvc[a]; (* an alias *) + +(* This is essentially a bracket of two operator expressions *) +scalarproductop[a_, b_] := vevwick[conj[a] ~ nc ~ zeroonvac[b]]; + +SetAttributes[{expvop, expvvc, braketop, braketvc}, Listable]; + +(* Helper function: this is a braket (creation operator rep) *) +braketop[a_, op_, b_] := Module[{ levo, desno }, + levo = conj[a] // Expand; + desno = op ~ nc ~ b // Expand; +(* desno = onvac @ desno; *) + vev[levo ~ nc ~ desno // Expand] + ]; + +(* Helper function: this is a braket (creation operator rep) *) +expvop[op_, a_] := braketop[a, op, a]; + +(* Helper function: this is a braket (occupation number rep) *) +braketvc[a_, op_, b_] := scalarproductvc[a, ap[op, b]]; + +(* An alias *) +braket[x___] := braketvc[x]; + +(* Helper function: this is a braket (occupation number rep) *) +expvvc[op_, a_] := braketvc[a, op, a]; + +(* Construct a matrix representation of operator 'a' in basis 'l'. + Note: 'l' must be expressed in occupation number representation. *) +matrixrepresentationvc[a_, l_List] := Module[{bl}, + bl = Map[ap[a, #]&, l]; + Outer[scalarproductvc, l, bl] +]; + +(* Construct a matrix representation of operator 'a' between two + different subspaces 'l1' and 'l2'. States 'l1' are bras, + while states 'l2' are kets. *) +matrixrepresentationvc[a_, l1_List, l2_List] := Module[{bl}, + bl = Map[ap[a, #]&, l2]; + Outer[scalarproductvc, l1, bl] +]; + +(* As above, but for operator expressions *) +matrixrepresentationop[a_, l_List] := Module[{bl}, + bl = Map[nc[a, #]&, l]; + Outer[scalarproductop, l, bl] +]; + +matrixrepresentationop[a_, l1_List, l2_List] := Module[{bl}, + bl = Map[nc[a, #]&, l2]; + Outer[scalarproductop, l1, bl] +]; + + +(************* ASCII <-> SNEG expression *************) + +(* a and b are delimiter strings, such as parenthesis. *) +snegToString[a_String, b_String] := StringJoin[a, b]; + +snegToString[a_String, b_String, i_] := StringJoin[a, ToString[i], b]; + +snegToString[a_String, b_String, l__] := + StringJoin[a, snegListToString[{l}], b]; + +snegListToString[l_List] := addcommas[ToString /@ l]; + +(* Concatenate a list of strings, separating them using commas *) +addcommas[l : {_String ..}] := + StringDrop[StringJoin[Map[# <> "," &, l]], -1]; + +(* Split a list of substrings separated by commas and transform them +to Mathematica expressions. *) +removecommas[s_String] := Map[ToExpression, StringSplit[s, ","]]; +removecommas[""] := {}; + +(* Replace operators, kets, bras with string expressions. Concatenate +the expression, taking into account the noncommutative multiplication +nc[]. *) +snegtoasciirules = { + op_[] /; operatorQ[op] :> ToString[op] <> "()", + op_[t_, i___] /; operatorQ[op] :> + ToString[op] <> If[t == CR, "+", ""] <> snegToString["(", ")", i], + ket[] :> "|>", + ket[i__] :> snegToString["|", ">", i], + bra[] :> "<|", + bra[i__] :> snegToString["<", "|", i], + HoldPattern[nc[l : _String ..]] :> StringJoin[l], + HoldPattern[nc[a___, s1_String, s2_String, b___]] :> + nc[a, StringJoin[s1, s2], b] +}; + +snegtoascii[expr_] := expr //. snegtoasciirules; + +(* Try to build an operator expression by parsing a string. *) + +buildop[op_String, type_] := With[{x = ToExpression[op]}, + x[type] /; (x =!= $Failed && operatorQ[x])]; + +buildop[op_String, type_, args_String] := + With[{x = ToExpression[op], arglist = removecommas[args]}, + x @@ Prepend[arglist, type] /; (x =!= $Failed && operatorQ[x])]; + +buildket[n_String] := ket @@ removecommas[n]; + +buildbra[n_String] := bra @@ removecommas[n]; + +snegstringparse[s_String] := StringCases[s, { + (op : Except[Characters["|<>()+"]] ..) ~~ + "+(" ~~ (args : Except[Characters[")"]] ...) ~~ ")" :> + buildop[op, CR, args], + (op : Except[Characters["|<>()+"]] ..) ~~ + "(" ~~ (args : Except[Characters[")"]] ...) ~~ ")" :> + buildop[op, AN, args], + "|" ~~ (ndx : Except[Characters["|<>()"]] ...) ~~ ">" :> + buildket[ndx], + "<" ~~ (ndx : Except[Characters["|<>()"]] ...) ~~ "|" :> + buildbra[ndx], + str__ :> str +}]; + +asciitosnegrules = { + s_String :> nc @@ snegstringparse[s] +}; + +asciitosneg[expr_] := expr //. asciitosnegrules; + +(************* (Q,Sz) BASIS **************) + +(* Make a basis of good quantum numbers Q and S_z. This is trivially done in +the occupation number representation vc[...]. *) +(* IMPORTANT: As of version 1.147, makebasis[] is called automatically. *) + +(* Spin-1/2 specific implementation! *) +qszbasisvc[l_List, 1/2] := Module[{n = Length[l], ndxs}, + makebasis[l]; + ndx2digs[ndx_] := IntegerDigits[ndx, 2, n]; + + qunum[up_, do_] := Module[{iup, ido}, + iup = Plus @@ ndx2digs[up]; + ido = Plus @@ ndx2digs[do]; + {iup+ido-n , 1/2 (iup-ido)} + ]; (* Q ans S_z *) + + ndx2vc[up_, do_] := vc @@ interleave[ndx2digs[up], ndx2digs[do]]; + + ndxs = Table[ {qunum[up, do], ndx2vc[up, do]}, + {up, 0, 2^n-1}, {do, 0, 2^n-1}]; + ndxs = Flatten[ndxs, 1]; + mergebasis[ndxs] +]; + +(* Generic implementation *) +qszbasisvc[l_List] := Module[{b, n, a, qmask, szmask}, + b = makebasis[l]; + n = Length[b]; + + a = Table[ IntegerDigits[i, 2, n], {i, 0, 2^n-1}]; + qmask = Table[1, {n}]; + szmask = b /. { + op_?fermionQ[CR, i___, SP_] /; spinof[op] == 1/2 :> + If[SP == UP, 1/2, -1/2], + op_?fermionQ[CR, i___, SP_] /; spinof[op] != 1/2 :> SP + }; + a = Map[{{qmask.# - n/2, szmask.#}, vc @@ #}&, a]; + a = mergebasis[a]; + a +]; + + +qszbasis[l_List, 1/2] := bzvc2bzop @ qszbasisvc[l, 1/2]; +qszbasis[l_List] := bzvc2bzop @ qszbasisvc[l]; + + +interleave[l__] := Flatten[ Transpose[{l}], 1]; + +(************* BASIS with no symmetries ***************) + +(* Make basis with no symmetries, i.e. no good quantum numbers. *) +nonebasisvc[l_List] := Module[{nr}, + makebasis[l]; + nr = Length[l]; + {{{}, Table[ vc @@ IntegerDigits[i, 2, 2nr], {i, 0, 4^nr-1}]}} +]; + +nonebasis[l_List] := bzvc2bzop @ nonebasisvc[l]; + +(********** BASIS with U(1) charge conservation symmetry **********) + +(* Make a basis with good quantum number Q *) + +qbasisvc[l_List] := Module[{nr, b}, + makebasis[l]; + nr = Length[l]; + b = Table[ vc @@ IntegerDigits[i, 2, 2nr], {i, 0, 4^nr-1}]; + b = Map[{{Count[#, 1, {1}] - nr}, {#}}&, b]; + b = mergebasis[b]; + b +]; + +qbasis[l_List] := bzvc2bzop @ qbasisvc[l]; + + +(************* (Q,S) BASIS ***************) + +(* Make a basis of good quantum numbers Q and S *) +(* Convention: always use the highest possible Sz, i.e. Sz == S ! *) + +qsbasisvcold[l_List] := bzop2bzvc[ qsbasis[l] ]; + + +qsbasisvc[l_List] := Module[{}, + makebasis[l]; + qsbasisvc[{}, l] +]; + +qsbasisvc[m_List, {}] := m; + +qsbasisvc[{}, l_List] := Module[{op1, op2}, + op1 = First[l] /. op_[j___] -> op[CR, j, UP]; + op2 = spinflip[op1]; + qsbasisvc[ + { + { {-1, 0}, {ap[1, VACUUM]} }, + { {0, 1/2}, {ap[op1, VACUUM]} }, + { {1, 0}, {ap[nc[op2, op1], VACUUM]} } + }, Drop[l, 1] + ] +]; + +qsbasisvc[m_List, l_List] := Module[{mmnew}, + mmnew = Simplify @ keymerge[ ADDvc[m, First[l]] ]; + qsbasisvc[mmnew, Rest[l]] +]; + +ADDvc[l_List, x:op_[___]] := + Map[ADDQSvc[#, x]&, l] ~ Flatten ~ 1; + +ADDQSvc[l_List, x:op_[___]] := + Map[ADDONEvc[#, x]&, Map[{First[l],#}&, l[[2]] ] ] ~ Flatten ~ 1; + +ADDONEvc[{{q_, s_}, a_}, op_[j___]] := Module[{op1, op2, res}, + op1 = op[CR, j, UP]; + op2 = spinflip[op1]; + res = { + { {q-1, s}, {a} }, + { {q+1, s}, {nc[op2, op1, a]} }, + { {q, s+1/2}, {nc[op1, a]} } + }; + If[s >= 1/2, Append[res, { + {q, s-1/2}, { Sqrt[2s/(2s+1)] nc[op2, a] - + Sqrt[1/(2s+1)] nc[op1, spindownvc[a]] } + }], res] +]; + +(* Change 0<->1 in the last index of an operator. This is just a faster +version of the more general invertspin[] function. *) +spinflip[op_[i_, j___, sigma_]] := op[i, j, 1-sigma]; + +(* We start with an empty list of (Q,S) states and the full list + of all single-particle operators. *) +(* makebasis is called automatically. While not necessary, such behave + is desirable for consistency with qszbasis[]. *) +qsbasis[l_List] := Module[{}, + makebasis[l]; + qsbasis[{}, l] +]; + +(* Last step: when no other operators need to be added, we're done! *) +qsbasis[m_List, {}] := m; + +(* First step: construct the three possible (Q,S) states from one + single-particle creation operator. *) +qsbasis[{}, l_List] := Module[{op1, op2}, + op1 = First[l] /. op_[j___] -> op[CR, j, UP]; + op2 = spinflip[op1]; + qsbasis[ + { + { {-1, 0}, {1} }, + { {0, 1/2}, {op1} }, + { {1, 0}, {nc[op2, op1]} } + }, Drop[l, 1] + ] +]; + +(* Intermediate steps: add one single-particle operator from list 'l' + and construct suitable combination of operators with a well-defined + (Q,S) quantum numbers. This step is reiterated until completion. *) +qsbasis[m_List, l_List] := Module[{mmnew}, + mmnew = Simplify @ keymerge[ ADD[m, First[l]] ]; + qsbasis[mmnew, Rest[l]] +]; + + +(********** BASIS with SU(2) spin invariance only **********) + +(* Approach: create a basis with (Q,S) quantum numbers, then +drop Q and finally merge subspaces with equal S quantum number. *) + +sbasisvc[l_List] := Module[{bvc}, + bvc = qsbasisvc[l]; + bvc = Map[{ {#[[1,2]]}, #[[2]] }&, bvc]; + bvc = mergebasis[bvc]; + bvc +]; + +sbasis[l_List] := bzvc2bzop @ sbasisvc[l]; + +(********** BASIS with U(1) spin invariance only **********) + +szbasisvc[l_List] := Module[{bvc}, + bvc = qszbasisvc[l]; + bvc = Map[{ {#[[1,2]]}, #[[2]] }&, bvc]; + bvc = mergebasis[bvc]; + bvc +]; + +szbasis[l_List] := bzvc2bzop @ szbasisvc[l]; + +(***** Miscelaneous functions used in basis construction *****) + +(* Returns 0 for operator expressions which annihilate vacuum *) +sneglinearoperatorFirst[zeroonvac]; +SetAttributes[zeroonvac, Listable]; + +zeroonvac[y:op_?fermionQ[___]] := + If[isannihilation[y], 0, y]; +zeroonvac[x:nc[a___, y:op_?fermionQ[i_, j___]]] := + If[isannihilation[y], 0, x]; + +zeroonvac[y:op_?bosonQ[AN,___]] := 0; +zeroonvac[x:nc[a___, y:op_?bosonQ[AN, j___]]] := 0; + +(* zeroonvac only works in the space of fermion operators *) +zeroonvac[k_ket] := k; +zeroonvac[nc[a___, k_ket]] := nc[zeroonvac[nc[a]], k]; + +(* If no rule was matched, we return the argument. *) +zeroonvac[a_] := a; + +nc[a___, x:op_?operatorQ[j___], VACUUM, b___] /; + (isannihilation[x]) := 0; +nc[a___, conj[VACUUM], x:op_?operatorQ[j___], b___] /; + (iscreation[x]) := 0; + +(* Vacuum state must be normalized to 1! *) +(* vev[nc[conj[VACUUM], VACUUM]] = 1; *) +nc[a___, conj[VACUUM], VACUUM, b___] := nc[a, b]; + +(* VACUUM placeholder corresponds to the vacuum state in the fermionic +space. Grassman numbers can therefore be moved to the left of the +conj[VACUUM] placeholder. *) +nc[a___, conj[VACUUM], z_?grassmanQ, b___] := + nc[a, z, conj[VACUUM], b]; + +If[PrettyOutput, + Format[VACUUM, StandardForm] := + StyleForm["\[LeftDoubleBracketingBar]\[Diameter]\[RightAngleBracket]", + FontColor -> Magenta]; + Format[conj[VACUUM], StandardForm] := + StyleForm["\[LeftAngleBracket]\[Diameter]\[RightDoubleBracketingBar]", + FontColor -> Magenta]; +]; + + +(* Smart spindown: determines which operators appear in the expressions, +builds a corresponding spin down operator and applies it to the expression. *) +SPINDOWN[a_] := nc[Plus @@ Map[spinminus, Union[Cases[a, +op_?operatorQ[i_, j___, sigma_] -> op[j], {0, Infinity} ]]], a]; + +(* Calculate a 'norm' for an operator expression, i.e. sqrt(). *) + +SetAttributes[{normop, normvc}, Listable]; + +normop[expr_] := Sqrt @ vev[ scalarproductop[expr, expr] ]; +normvc[expr_] := Sqrt @ vev[ scalarproductvc[expr, expr] ]; + +(* applytocr[] applys an operator to an expression, drops parts which +annihilate the vacuum and "normalizes" the reminder *) +applytocr[op_, expr_] := Module[{b}, + b = zeroonvac[nc[op, expr]]; + b / normop[b] +]; + +(* Smarter spindown: drops parts which annihilate vacuum and "normalizes" +the remainder *) +spindown::Usage = "spindown[v] effectively applies the spin-lowering +operator to the state v in the creation operator representation."; +spindown[a_] := Module[{b}, + b = zeroonvac[SPINDOWN[a]]; + b / normop[b] +]; + + +(* spindown function for states in the occupation number representation *) +spindownvc::Usage = "spindownvc[v] effectively applies the spin-lowering +operator to the state v in the occupation number representation."; + +sneglinearoperator[spindownvc1]; +spindownvc1[a_vc] := Module[{b, fn}, + b = List @@ Partition[a,2]; + fn[vc[1,0], {i_}] := vc @@ (Join @@ MapAt[vc[0,1]&, b, i]); + fn[_, {_}] := 0; + Total[MapIndexed[fn, b]] +]; + +spindownvc[a_] := spindownvc1[a] // (If[#=!=0, #/normvc[#], 0])&; + + +keyselect[l_, key_] := Select[l, First[#] == key &]; + +keymerge[l_] := Module[{keys, sublists}, + keys = Union[ l[[All,1]] ]; + sublists = Map[keyselect[l, #]&, keys]; + Table[{keys[[i]], sublists[[i,All,2]] // Flatten}, {i, Length[keys]}] +]; + +(* Add one (doublet) operator to a (Q,S) basis. This is used + by qsbasis[] to generate the entire list of basis states. *) +ADD[l_List, x:op_[j___]] := + Map[ADDQS[#, x]&, l] ~ Flatten ~ 1; + +ADDQS[l_List, x:op_[j___]] := + Map[ADDONE[#, x]&, Map[{First[l],#}&, l[[2]] ] ] ~ Flatten ~ 1; + +ADDONE[{{q_, s_}, a_}, op_[j___]] := Module[{op1, op2, res}, + op1 = op[CR, j, UP]; + op2 = spinflip[op1]; + res = { + { {q-1, s}, {a} }, + { {q+1, s}, {nc[op2, op1, a]} }, + { {q, s+1/2}, {nc[op1, a]} } + }; + If[s >= 1/2, Append[res, { + {q, s-1/2}, { Sqrt[2s/(2s+1)] nc[op2, a] - + Sqrt[1/(2s+1)] nc[op1, spindown[a] ] } + }], res] (* cf. Krishna-Murthy et al. *) +]; + +(**** Projection operators ****) + +(* Indexing is the same as in T. Rejec's operachain.cpp code for + evaluation of operator expressions *) + +(* Short-cuts *) +ClearAttributes[{PROJ0, PROJUP, PROJDO, PROJ2, PROJ1, PROJ02}, Protected]; +PROJ0 = 1; +PROJUP = 2; +PROJDO = 3; +PROJ2 = 4; +PROJ1 = 5; +PROJ02 = 6; +SetAttributes[{PROJ0, PROJUP, PROJDO, PROJ2, PROJ1, PROJ02}, Protected]; + +SetAttributes[{projector, projector0, projectorUP, projectorDO, +projector1, projector2, projector02}, Listable]; + +(* Subspace with 0 particles *) +projector[op_[j___], 1] := nc[ 1-number[op[j], UP], 1-number[op[j], DO] ]; +projector0[op_[j___]] := projector[op[j], PROJ0]; + +(* Subspace with spin up particle *) +projector[op_[j___], 2] := nc[ number[op[j], UP], 1-number[op[j], DO] ]; +projectorUP[op_[j___]] := projector[op[j], PROJUP]; + +(* Subspace with spin down particle *) +projector[op_[j___], 3] := nc[ number[op[j], DO], 1-number[op[j], UP] ]; +projectorDO[op_[j___]] := projector[op[j], PROJDO]; + +(* Subspace with two particles *) +projector[op_[j___], 4] := nc[ number[op[j], UP], number[op[j], DO] ]; +projector2[op_[j___]] := projector[op[j], PROJ2]; + +(* Subspace with one particle of either spin *) +projector[op_[j___], 5] := projector[op[j], 2] + projector[op[j], 3]; +projector1[op_[j___]] := projector[op[j], PROJ1]; + +(* Subspace with occupancy 0 or 2 *) +projector[op_[j___], PROJ02] := + projector[op[j], PROJ0] + projector[op[j], PROJ2]; +projector02[op_[j___]] := projector[op[j], PROJ02]; + +projector[op_?AtomQ, n___] := projector[op[], n]; +projector0[op_?AtomQ] := projector0[op[]]; +projectorUP[op_?AtomQ] := projectorUP[op[]]; +projectorDO[op_?AtomQ] := projectorDO[op[]]; +projector1[op_?AtomQ] := projector1[op[]]; +projector2[op_?AtomQ] := projector2[op[]]; +projector02[op_?AtomQ] := projector02[op[]]; + + +(**** Operations on basis sets ****) +bzvc2bzop[bvc_] := Map[{First[#], vc2ops[#[[2]]]}&, bvc]; +(* If the vacuum state vac is not provided, vacuum[] is used. *) +bzop2bzvc[bz_, vak_] := Map[{First[#], ap[#[[2]], vak]}&, bz]; +bzop2bzvc[bz_] := bzop2bzvc[bz, vacuum[] ]; + +(* Decompose a vector in a given basis *) +decomposeop[vec_, basis_] := Map[scalarproductop[#, vec]&, basis]; +decomposevc[vec_, basis_] := Map[scalarproductvc[#, vec]&, basis]; + + +(**** Orthogonalization ****) + +orthogop[vecs_?VectorQ, basis_?VectorQ] := Module[{m, l}, + m = Map[ decomposeop[#, basis]&, vecs ]; + m = Simplify[m]; (* See below! *) + l = snegorthog[m]; + l = Simplify[l]; + l = Select[l, Norm[#] != 0 &]; (* Drop zero vectors *) + Expand @ Map[(# . basis)&, l] +]; + +(* TO DO: bottleneck! Can we do it in a more efficient manner? *) +orthogvc[vecs_?VectorQ, basis_?VectorQ] := Module[{m, l}, + m = Map[ decomposevc[#, basis]&, vecs ]; + + (* Occasionally, when dealing with very complex sets of basis vectors, + decomposevc[] returns large non-simplified expressions. We need to + Simplify them before performing the orthogonalisation. *) + + m = Simplify[m]; + l = snegorthog[m]; + l = Simplify[l]; + l = Select[l, Norm[#] != 0 &]; (* Drop zero vectors *) + Expand @ Map[(# . basis)&, l] +]; + +(*** Projecting out a state ***) + +(* Project out 'state' from a space of 'vecs'. Is is assumed that + both 'state' and 'vecs' are normalized to 1. + The resulting vectors are reorthogonalized. *) +remove1state[vecs_, state_] := Module[{vecs2}, + If[vecs != {}, + vecs2 = Map[# - scalarproductvc[state, #] state &, vecs]; + vecs2 = orthogvc[vecs2, vecs] // Expand, + vecs2 = {} + ]; + vecs2 +]; + +(* Project out 'states' from a space of 'vecs'. *) +(* XXX: preformance bottleneck for transform1[]!!! *) +removestates[vecs_, states_] := Module[{}, + Fold[remove1state, vecs, states] +]; + +(** +Transformation from basis with well defined (Q,S) to a basis with well +defined (I,S) quantum numbers. + +Strategy: take states with high Q (these are maximum weight states with +respect to the isospin SU(2) group), apply Tminus to generate other states +of the multiplets, and subtract those states from the lower Q subspaces. +Repeat until no low Q states remain. At this point, the basis consists only +of maximum-weight representative states, which is the (I,S) basis that we +wanted to construct. + +NOTE: this function may also be applied to a (Q,Sz) basis to generate a +(I,Sz) basis. In fact, it is applicable to all (Q,...) -> (I,...) +transformations. +**) + +transformQStoIS[baza_, TM_] := transformQStoIS2[TM, Reverse[baza], {}]; + +transformQStoIS::Tminus = "Operator Tminus is not defined."; +transformQStoIS[baza_] := Module[{}, + If[!ValueQ[Tminus], Message[transformQStoIS::Tminus]; Return[$Failed]]; + transformQStoIS2[Tminus, Reverse[baza], {}] +]; + +transformQStoIS2[TM_, {}, novabaza_] := dropemptysubspaces[novabaza]; +transformQStoIS2[TM_, baza_, novabaza_] := Module[{nova2, ISv, ostali}, + {ISv, ostali} = transform1[baza, TM]; + nova2 = Append[novabaza, ISv]; + transformQStoIS2[TM, ostali, nova2] +]; + +(* baza is a basis List in the format as returned by function qsbasis[] *) +(* Tminus is the isospin lowering operator that we apply to the vectors in + the subspace that we're subtracting from other subspaces *) +transform1[baza_, Tminus_] := + Module[{prvi, ostali, subspace, states, statesplus, pos, ISvec, i, + nrqn, subvector, isubspace}, + prvi = First[baza]; (* The subspace that we are subtracting *) + ostali = Rest[baza]; (* The subspaces that we are subtracting from *) + + (* subspace=(Q, ...) indexes. Note that no assumptions is made about + the indexes past the first one. *) + subspace = First[prvi]; + nrqn = Length[subspace]; (* Total number of quantum numbers *) + (* (1,0,...) vector, used in subtraction of quantum numbers later on *) + subvector = Table[If[i==1, 1, 0], {i, nrqn}]; + + (* Q -> I change of the first quantum number *) + isubspace = subspace; isubspace[[1]] = (1/2) isubspace[[1]]; + + states = prvi[[2]]; (* Vectors *) + + (* Vectors to which we are consecutively applying Tminus *) + statesplus = states; + + For[i = 1, True, i++, + statesplus = Map[ap[Tminus, #] &, statesplus]; + statesplus = Expand @ statesplus; (* WAS: Simplify *) + statesplus = DeleteCases[statesplus, 0]; (* Drop nul vectors *) + If[statesplus == {}, + Break[]; (* We're done *) + ]; + + (* Normalize to 1 *) + statesplus = Map[#/normvc[#] &, statesplus]; + statesplus = Expand @ statesplus; (* WAS: Simplify *) + + (* Find relevant subspaces from which we can project out *) + pos = Position[ostali, x_ /; First[x] == subspace - (2i)subvector, {1}, + Heads -> False]; + + If[pos == {}, + Break[], + (* Performing the projecting *) + ostali = MapAt[{First[#], removestates[#[[2]], statesplus]} &, + ostali, pos]; + ]; + ]; + + (* Return {(I, ...) pair, vectors} and the matrix of remaining vectors *) + ISvec = {isubspace, states}; + {ISvec, ostali} +]; + +(*** BRAKET ***) + +(* Dirac's bra-ket notation for arbitrary states. This is useful to represent +non-fermionic degrees of freedom, such as bosons, spins, rotators, etc. *) + +If[PrettyOutput, + (* Enclose bras and kets in suitable angle brackets and bars (orange!), + replace Null elements by small circles. *) + PrettyOutputBraKetRule = { Null -> "\[SmallCircle]" }; + + Format[bra[i___], StandardForm] := + DisplayForm @ RowBox[Flatten[{ + StyleForm["\[LeftAngleBracket]", FontColor -> Orange], + If[{i} =!= {}, RowBox[{i} /. PrettyOutputBraKetRule], {}], + StyleForm["\[RightBracketingBar]", FontColor -> Orange] + }, 1]]; + Format[ket[i___], StandardForm] := + DisplayForm @ RowBox[Flatten[{ + StyleForm["\[LeftBracketingBar]", FontColor -> Orange], + If[{i} =!= {}, RowBox[{i} /. PrettyOutputBraKetRule], {}], + StyleForm["\[RightAngleBracket]", FontColor -> Orange] + }, 1]]; + + fixparens[bra]; + fixparens[ket]; +]; + +(* Convention: always commute bras and kets to the right of the +fermionic (!) operators *) + +nc[a___, x1:bra[b___], x2:op_?fermionQ[i___], d___] := nc[a, x2, x1, d]; +nc[a___, x1:ket[k___], x2:op_?fermionQ[i___], d___] := nc[a, x2, x1, d]; + +(* Returns a list with True at non-Null positions *) +nonnullpattern[l_List] := nonnullQ /@ l; +nonnullQ[expr_] := (expr =!= Null); + +(* Returns True if there are no conflicting non-Null elements in two lists. This +is used in concatenating kets and bras. *) +compatiblepattern[l1_List, l2_List] := (Length[l1] == Length[l2]) && (Nor @@ + Thread[And[nonnullpattern[l1], nonnullpattern[l2]]]); + +(* Returns True if non-Null elements come in pairs. This is used in calculating +scalar product, i.e. expectation value. *) +pairpattern[l1_List, l2_List] := (Length[l1] == Length[l2]) && !(Or @@ + Thread[Xor[nonnullpattern[l1], nonnullpattern[l2]]]); + +(* Return position that do not contain the Null element *) +nonnullpos[l_List] := Position[l, _?nonnullQ, {1, Infinity}, Heads->False]; + +(* Default rule for the expectation values: Kronecker delta over all non-Null indexes *) +braketrule[bra[b___], ket[k___]] /; pairpattern[{b}, {k}] := + Module[{pos, l1, l2}, + pos = nonnullpos[{b}]; + l1 = Extract[{b}, pos]; + l2 = Extract[{k}, pos]; + Inner[KroneckerDelta, l1, l2, Times] + ]; + +(* When bra and ket are juxtaposed, we can evaluate the scalar +product. We only do this when the bra and ket are compatible in the +sense that both have the same pair pattern: quantum numbers at +equivalent position are defined either in both bra and ket or in +neither. *) + +nc[a___, x1:bra[b___], x2:ket[k___], d___] /; pairpattern[{b}, {k}] := + braketrule[x1, x2] nc[a, d] + +(* bkcombine[]: combine two lists in such a manner that Null elements +are treated as placeholders that are filled only if the corresponding +element in the other list is non-Null. Used in concatenation of bras +and kets.*) +(* TO DO: This could be made to work even for multi-dimensional +lists. *) + +bkcombine::incompat = "Incompatible pattern of non-Null elements"; +bkcombine[l1_List, l2_List] := Module[{pos}, + pos = nonnullpos[l2]; + If[pos =!= {}, ReplacePart[l1, l2, pos, pos], l1] +] /; compatiblepattern[l1, l2]; + +bkcombine[l1_List, l2_List] := (Message[bkcombine::incompat]; + error[bkcombine, l1,l2]) /; !compatiblepattern[l1, l2]; + +(* Direct product: concatenates all indexes, if possible! *) +nc[a___, bra[i1___], bra[i2___], b___] /; + compatiblepattern[{i1}, {i2}] := nc[a, bra @@ bkcombine[{i1}, {i2}], b]; +nc[a___, ket[i1___], ket[i2___], b___] /; + compatiblepattern[{i1}, {i2}] := nc[a, ket @@ bkcombine[{i1}, {i2}], b]; + +(* Trick: bra and ket from orthogonal Hilbert spaces: we may commute them! *) +nc[a___, x1:bra[i1___], x2:ket[i2___], b___] /; + compatiblepattern[{i1}, {i2}] := nc[a, x2, x1, b]; + +(* Conjugation *) +conj[ket[k___]] := bra[k]; +conj[bra[b___]] := ket[b]; + +(* Build an object with head 'head', and Null elements at positions where +'positionarray' contains zeros, and 'value' at positions where +'positionarray' contains ones. *) +snegbuild[head_, value_, positionarray_] := + head @@ (positionarray /. {0 -> Null, 1 -> value}); +buildket[qn_, positionarray_] := snegbuild[ket, qn, positionarray]; +buildbra[qn_, positionarray_] := snegbuild[bra, qn, positionarray]; + +(* Spin states and spin operators *) +spinket[s_?halfintegerQ, positionarray_:{1}] := + buildket[#, positionarray]& /@ Range[s, -s, -1]; +spinbra[s_?halfintegerQ, positionarray_:{1}] := + buildbra[#, positionarray]& /@ Range[s, -s, -1]; + +spinketbraZ[s_?halfintegerQ, positionarray_:{1}] := + VMV[spinket[s, positionarray], spinmatrixZ[s], spinbra[s, positionarray]]; +spinketbraP[s_?halfintegerQ, positionarray_:{1}] := + VMV[spinket[s, positionarray], spinmatrixP[s], spinbra[s, positionarray]]; +spinketbraM[s_?halfintegerQ, positionarray_:{1}] := + VMV[spinket[s, positionarray], spinmatrixM[s], spinbra[s, positionarray]]; +spinketbraX[s_?halfintegerQ, positionarray_:{1}] := + VMV[spinket[s, positionarray], spinmatrixX[s], spinbra[s, positionarray]]; +spinketbraY[s_?halfintegerQ, positionarray_:{1}] := + VMV[spinket[s, positionarray], spinmatrixY[s], spinbra[s, positionarray]]; + +(* Identity operator *) +spinketbraI[s_?halfintegerQ, positionarray_:{1}] := + VMV[spinket[s, positionarray], IdentityMatrix[2s+1], + spinbra[s, positionarray]]; + +(* Convenience definitions *) +spinketbraOPS[s_?halfintegerQ, positionarray_: {1}] := + (#[s, positionarray]) & /@ + {spinketbraX, spinketbraY, spinketbraZ, spinketbraP, spinketbraM}; +spinketbraXYZ[s_?halfintegerQ, positionarray_: {1}] := + (#[s, positionarray]) & /@ + {spinketbraX, spinketbraY, spinketbraZ}; + +(* Spin-spin exchange coupling *) +spinketbraspinspin[{s1_?halfintegerQ, s2_?halfintegerQ}] := + Module[{sx1, sy1, sz1, sx2, sy2, sz2}, + {sx1, sy1, sz1} = spinketbraXYZ[s1, {1, 0}]; + {sx2, sy2, sz2} = spinketbraXYZ[s2, {0, 1}]; + Expand[sx1~nc~sx2 + sy1~nc~sy2 + sz1~nc~sz2] + ]; + +(* Generates a basis for spin operators. The key is the value of + the Z-component of spin, S_z. *) +spinbasis[s_?halfintegerQ, positionarray_:{1}] := + {{#}, {buildket[#, positionarray]}}& /@ Range[s, -s, -1]; + +(* Generates a direct-product basis for two spin operators. The key + is the Z-component of the TOTAL spin. *) +spinbasis[{a_?halfintegerQ, b_?halfintegerQ}] := + basistensorproduct[spinbasis[a, {1, 0}], + spinbasis[b, {0, 1}], + Function[{sz1, sz2}, sz1+sz2] ]; + +spinbasis[{a_?halfintegerQ, b__}] := + basistensorproduct[spinbasis[a, {1, 0}], + spinbasis[b, {0, 1}], + Function[{sz1, sz2}, sz1+sz2] ]; + +(* Apply rule to the basis states (vektor form) in 'baza'. +Used by transformbasis[]. *) +transformfunc[baza_, rule_] := Module[{vecs}, + vecs = Map[Simplify[ ap[rule, #] ]&, baza]; + orthogvc[vecs, baza] +]; + +(* Apply rule to an entire basis in the form of a list of {{Q. numbers}, +{basis states}}... *) +transformbasis[bvc_List, rule_] := Module[{bvcnew}, + bvcnew = Map[{First[#], transformfunc[ #[[2]], rule ]}&, bvc]; + dropemptysubspaces[bvcnew] +]; + +(* Apply several rules to each basis element. *) +(* XXX: Untested! *) +transformbasismulti[bvc_List, rules_List] := Module[{bvcnew}, + bvcnew = Map[transformbasis[bvc, #]&, rules]; + bvcnew = Flatten[bvcnew, 1]; + mergebasis[bvcnew] +]; + +(* Merge the basis: gather basis states with the same "key". *) +mergebasis[bvc_List] := Module[{qn}, + qn = Union @ bvc[[All,1]]; (* Get all different quantum numbers. *) + Map[{#, collectvalues[#, bvc]}&, qn] +]; + +(* Make a spinless basis by projecting out all states with spin-down +electrons. *) +spinlessbasis[bzin_List] := Module[{bz}, + bz = bzin //. { op_?fermionQ[CR, i___, DO] -> 0 }; + bz = Map[{First[#], DeleteCases[Last[#], 0]}&, bz]; + bz = mergebasis[bz]; + bz = DeleteCases[bz, {{___},{}}] +]; + +(* Get all basis vectors with the same quantum numbers (key). *) +collectvalues[key_, bvc_List] := Flatten @ + (Select[bvc, (First[#] === key)&] [[All,2]]); + +(* Drop subspaces with no states from a given basis bvc *) +dropemptysubspaces[bvc_List] := Select[bvc, #[[2]] =!= {} &]; + +(* Apply a function to an entire basis. *) +applybasis[bvc_?bzQ, fn_] := Map[{First[#], fn /@ Last[#]}&, bvc]; + +(* Note: spin basis in S representation can be conveniently +constructed using qsbasis[] which includes a pseudofermion that +corresponds to the spin, and then projecting out zero and double +occupancy states, and finally replacing the remaining single particle +pseudofermion operators with spin kets. One should be careful about +the anticommutation operators, though. The pseudospin operator must be +moved to the rightmost position of the operator string by overriding +the conventional (anti)commutation rules. Of course, this only works +for spin 1/2. *) + +QSpseudospin[bvc_List, x:op_?operatorQ[]] := Module[{vecs}, + vecs = transformbasis[bvc, projector[x, PROJ1]]; + vecs = applybasis[vecs, vc2ops]; + vecs = vecs /. op[j___] :> hold[op[j]]; + vecs = vecs //. nc[a___, x1:hold[op[j1___]], x2:b_?operatorQ[j2___], c___] :> + -nc[a, x2, x1, c]; + vecs = vecs /. hold[op[i_, sp_, j___]] :> ket[If[sp == UP, 1/2, -1/2]] +]; + +(* basistensorproduct combines two basis sets bz1 and bz2 by forming +tensor products. Function qnfn defines how to combine the lists of +quantum numbers. *) +basistensorproduct[bz1_List, bz2_List, qnfn_:List] := Module[{fn}, + fn[{qn1_, vec1_}, {qn2_, vec2_}] := + {qnfn[qn1, qn2], Flatten @ outer[vec1, vec2]}; + mergebasis[Flatten[Outer[fn, bz1, bz2, 1], 1]] +]; + +(* Create a new basis by forming products of basis states from each of the + basis given as arguments. qnmergefnc specifies how the lists of quantum numbers + are to be merged to give the quantum number corresponding to the subspaces + of the new basis. *) +(* XXX: See also basistensorproduct[] - essentially the same! *) +directproductbasis[basis1_List, basis2_List, qnmergefnc_:Join] := + Flatten[Outer[directproductsubspace[##, qnmergefnc] &, basis1, basis2, 1], 1]; + +(* Perform the merging with the pair of basis subspaces. *) +directproductsubspace[{qn1_List, states1_List}, {qn2_List, states2_List}, + qnmergefnc_] := {qnmergefnc[qn1, qn2], Flatten[outer[states1, states2], 1]}; + + +(* Add a spin-ket for spin 'SP' to a (Q,Sz) basis. *) +QSZaddspinket[bvc_List, SP_] := Module[{spinbz}, + spinbz = spinbasis[SP]; + basistensorproduct[spinbz, bvc, Function[{qn1, qn2}, + {First[qn2], Last[qn2] + First[qn1]}]] +]; + +(* Application of an operator to a occupation number representation vector *) +ap[x:ket[___], vc[]] := vc[x]; +ap[x:ket[___], vc[j__]] := vc[j, x] /; Head[Last[{j}]] =!= ket; +ap[ket[i___], vc[j___, ket[k___]]] := vc[j, ket @@ bkcombine[{i}, {k}]] /; + compatiblepattern[{i}, {k}]; +ap[a___, x1:bra[i___], vc[j___, x2:ket[k___]]] := ap[a, vc[j]] braketrule[x1, x2] /; + pairpattern[{i}, {k}]; + +(* Direct product of a fermionic basis and a (single) phonon basis +with up to Nph phonons. *) + +transformtoPH[bvc_List, Nph_Integer] := Module[{vecs, fn, phbasis}, + phbasis = phononbasis[Nph]; + fn[x_] = nc[x, #]&; + vecs = Map[applybasis[bvc, fn[#]]&, phbasis]; + mergebasis @ Flatten[vecs, 1] +]; + +phononbasis[Nph_] := Table[ket[i], {i, 0, Nph}]; + +(* Phonon hamiltonian construction functions *) +phononnumber[Nph_Integer] := Sum[i nc[ket[i], bra[i]], {i, 0, Nph}]; +phononplus[Nph_Integer] := Sum[ Sqrt[i+1] nc[ket[i+1], bra[i]], {i, 0, Nph-1}]; +phononminus[Nph_Integer] := Sum[ Sqrt[i] nc[ket[i-1], bra[i]], {i, 1, Nph}]; +phononx[Nph_Integer] := phononplus[Nph] + phononminus[Nph]; + +(*** LEFT/RIGHT SYMMETRY SUPPORT ***) + +(* Given l, the list of operators, lrmap[] returns the transformation rules +that perform the left/right transformation. Since mbfunc[] has attribute +Listable, l can be a nested list. In this case, each sublist denotes sites +which are equivalent under the mirror symmetry. For example: {a[], {b[], +c[]}, d[]}: in this case, a maps to d and vice versa, but c and d are +invariant. Or for {{a[],b[]},{c[],d[]}}: a maps to c, b maps to d, and vice +versa. *) + +lrmap[l_List] := Module[{l1, l2}, + l1 = Flatten @ (mbfunc /@ l); + l2 = Flatten @ (mbfunc /@ Reverse[l]); + Thread[l1 -> l2] +]; + +(* Perform an operation o on each subspace vectors *) + +mapbasis::compat = "Incompatible subspaces"; +mapbasis[bz1_, bz2_, o_:(1/2 Plus[##]&)] := Module[ + {subs1, vecs1, subs2, vecs2, vecs}, + subs1 = bz1[[All, 1]]; + vecs1 = bz1[[All, 2]]; + subs2 = bz2[[All, 1]]; + vecs2 = bz2[[All, 2]]; + If[subs1 =!= subs2, + Message[mapbasis::compat]; + Return[]; + ]; + + vecs = Table[o[vecs1[[i]], vecs2[[i]]], {i, Length[vecs1]}]; + + Transpose[{subs1, vecs}] +]; + +(* Orthogonalize bz1 with respect to original basis bz *) +orthogbasisop[bz1_, bz_] := mapbasis[bz1, bz, orthogop]; +orthogbasisvc[bvc1_, bvc_] := mapbasis[bvc1, bvc, orthogvc]; + +(* Append a quantum number 'qn' to the good quantum numbers of each subspace +in basis bz *) +appendquantumnumber[bz_, qn_] := Map[{ Append[First[#], qn], #[[2]] }&, bz]; + +(* Putting it all together: transform a given basis to an equivalent basis +which takes into account the l/r symmetry, i.e. parity. l is the list of +consecutive operators that are fed to lrmap[] function. Additional +transformation rules can be specified as an optional parameter +'snextrarule'.*) +(* The input basis 'bz' must be in the operator form. *) + +transformtoLR[bz_, l_, snextrarule_:{}] := Module[{bz2, bzsim, bzasim}, + bz2 = bz /. Join[lrmap[l], snextrarule]; + bzsim = mapbasis[bz, bz2, 1/2 Plus[##]& ]; + bzasim = mapbasis[bz, bz2, 1/2 Subtract[##]& ]; + {bzsim, bzasim} = Map[dropemptysubspaces@orthogbasisop[#, bz]&, {bzsim, bzasim}]; + bzsim = appendquantumnumber[bzsim, 1]; + bzasim = appendquantumnumber[bzasim, -1]; + Sort @ Join[bzsim, bzasim] +]; + +transformtoLRvc[bz_, l_, vak_, snextrarule_:{}] := Module[ + {bvc, bz2, bzsim, bzasim, bvcsim, bvcasim}, + bvc = bzop2bzvc[bz, vak]; + bz2 = bz /. Join[lrmap[l], snextrarule]; + bzsim = mapbasis[bz, bz2, 1/2 Plus[##]& ]; + bzasim = mapbasis[bz, bz2, 1/2 Subtract[##]& ]; + bvcsim = bzop2bzvc[bzsim, vak]; + bvcasim = bzop2bzvc[bzasim, vak]; + {bvcsim, bvcasim} = + Map[dropemptysubspaces@orthogbasisvc[#, bvc]&, {bvcsim, bvcasim}]; + bvcsim = appendquantumnumber[bvcsim, 1]; + bvcasim = appendquantumnumber[bvcasim, -1]; + Sort @ Join[bvcsim, bvcasim] +]; + + +(** Transformations of basis sets **) + +(* Produce a vector with 1 on i-th site in a length-n vector *) +enain[i_, n_] := Table[If[j == i, 1, 0], {j, n}]; + +(* Given a set of operators "oldbasis", a set of operators "newbasis" and +a set of operator expressions for newbasis, produce the transformation +rules from the oldbasis to the new basis! *) +snegold2newrules[oldbasis_, newbasis_, rules_] := Module[{n, ob, mat, x}, + n = Length[oldbasis]; + If[n != Length[newbasis], + Print["Mismatching lengths."]; + Return[]; + ]; + ob = Table[oldbasis[[i]] -> enain[i, n], {i, Length[oldbasis]}]; + mat = rules /. ob; + x = Norm[mat . Transpose[Conjugate[mat]] - IdentityMatrix[n]]; + If[x != 0, + Print["The transformation rule is not unitary."]; + Return[]; + ]; + mat = Inverse[mat]; + Thread[Rule[oldbasis, mat . newbasis]] +]; + + +(*** WICK'S THEOREM ***) + +(* Produce all possible pairs with elements taken from two lists *) +(* OPTIMIZATION: cache the result! *) +allpairs[m_List, n_List] := allpairs[m,n] = Flatten[Outer[List, m, n], 1]; + +(* Return a list of positions of CR/AN opperators *) +indecesCR[HoldPattern[nc[l__]]] := Flatten @ Position[{l}, op_[CR, j___], {1}]; +indecesAN[HoldPattern[nc[l__]]] := Flatten @ Position[{l}, op_[AN, j___], {1}]; + +(* Return the maximal number of pairs that can be contracted. *) +maxcran[HoldPattern[x : nc[l__]]] := + Min[Length @ indecesCR[x], Length @ indecesAN[x]]; + +(* Returns True if there are no overlapping elements in the list of pairs *) +disjointQ[{l__}] := + Length[{l}] == 1 || Length[Union[l]] == (Plus @@ Map[Length, {l}]); + +(* Conflicting lists of pairs are zeroed *) +zeroconflict[z_List] := Map[If[disjointQ[#], #, 0] &, z]; + +(* Drop 0 elements from a List *) +dropzeros[l_List] := Delete[l, Position[l, 0, {1}]]; + +(* Drop lists of pairs that have overlapping indecs *) +dropconflicts[l_List] := dropzeros @ zeroconflict[l]; + +(* Generate lists of nr pairs to be contracted *) +(* OPTIMIZATION: cache the result! *) +contractedpairs[nr_, l_List] := contractedpairs[nr,l] = + dropconflicts @ KSubsets[l, nr]; + +(* Contraction of two operators according to the standard definition: +substract the normal ordered expression. *) +contraction[x1 : op1_[t1_, i1___], x2 : op2_[t2_, i2___]] /; + (ordering[op] =!= NONE) := + nc[x1, x2] - normalorder[x1 ~ nc ~ x2]; + +(* Contract one pair. ctr is an order index that replaces the contracted +operators in the operator string, factor is the prefactor of the string and +ops is the string that we are currently applying the contraction to. *) + +contractone[{ctr_, factor_, ops_List}, {i1_, i2_}] := + Module[{l = ops, c}, + (* Preserve order! *) + If[i1 < i2, + c = contraction[l[[i1]], l[[i2]]], + c = contraction[l[[i2]], l[[i1]]] + ]; + l[[i1]] = ctr; + l[[i2]] = ctr; + {ctr + 1, c*factor, l} + ]; + +(* Return the remaining operator string. The integer number indexed +placeholders are commuted to the beginning of the list and ordered: +this gives the correct sign of the remaining operator string. Sign is +flipped whenever two neighboring integers are transposed, and when an +integer and a neighboring operator are transposed. Note: List2 is a +non-Listable List. *) + +remainder[{ops__}] := Module[{ops2}, + ops2 = List2[ops] //. { + List2[a___, nr1_Integer, nr2_Integer, b___] :> + -List2[a, nr2, nr1, b] /; nr2 < nr1, + List2[a___, x : c_[i__], nr_Integer, b___] :> + -List2[a, nr, x, b] /; operatorQ[c] + }; + ops2 = ops2 /. {x : List2[i__] :> DeleteCases[x, _Integer, {1}]}; + ops2 = ops2 /. {List2[i___] :> dd[nc[i]]} +]; + +(* Contract "pairs" in the operator string "ops" *) + +contract[ops_List, pairs_List] := Module[{ctr, factor, indexes}, + {ctr, factor, indexes} = Fold[contractone, {1, 1, ops}, pairs]; + remainder[indexes] * factor +]; + +(* Do the Wick! *) + +wickorder[0, HoldPattern[x : nc[l__]]] := dd[x]; + +wickorder[nr_ /; nr > 0, HoldPattern[x : nc[l__]]] := With[ + { + cp = contractedpairs[nr, allpairs[indecesCR[x], indecesAN[x]]] + }, + Plus @@ Map[contract[{l}, #] &, cp] +]; + +sneglinearoperator[wick]; +SetAttributes[wick, Listable]; + +(* Write an operator string using normal ordered strings and contractions *) +wick[HoldPattern[x : nc[l__]]] := Module[{i, m}, + m = maxcran[x]; + Sum[wickorder[i, x], {i, 0, m}] +]; +wick[] := 1; + +(* wick for formal sums *) +wick[sum[a_, b_]] := sum[wick[a], b]; + +(* VEV calculation using the Wick theorem *) +(* Idea: 1) the expression is rewritten in terms of contractions and +normal ordered terms. 2) vev is taken. 3) All terms but the last one +(with nr=maxcran[x]) drop out since vev[dd[...]]=0. In other words, +we do the same calculation as wick[] does, but calling wickorder[] +for a single argument. *) + +vevwick[HoldPattern[x : nc[l__]]] := wickorder[maxcran[x], x] /. {dd[_]->0}; +vevwick[] := 1; + +(* vevwick2[] is a more efficient implementation of vevwick[]: it +immediately drops terms which are zero. Note that here we made certain +assumptions about the vacuum state and the normal ordered products, +thus this version of vevwick is not fully general! It is much faster, +though. Added 4.3.2010 *) + +sneglinearoperatorFirst[vevwick2]; +vevwick2[HoldPattern[x : nc[l__]]] := Module[{ic, ia, ndc, nda}, + ic = indecesCR[x]; + ia = indecesAN[x]; + If[Length[ic] != Length[ia], Return[0]]; + (* Bug trap *) + ndc = Check[Extract[x, Map[{#} &, ic]] /. op_[CR, i___] -> {op, i}, + Null]; + nda = Check[Extract[x, Map[{#} &, ia]] /. op_[AN, i___] -> {op, i}, + Null]; + If[ndc === Null || nda === Null, Return[Problem[{l}]] ]; + If[Sort[ndc] =!= Sort[nda], Return[0]]; + (* The term contributes. Perform the calculation! *) + vevwick[x] +]; + + +(* Generate an (I,S) basis, for example: quickISObasis[{a[],b[],c[]}]. +NOTE: nnop[] will be used, if defined, otherwise all signs are taken to be +equal (this is likely not what you want!). *) + +quickbasis[basisops__, fnc_] := Module[{Tminus, nnop2}, + nnop2[i_] := If[ValueQ @ nnop[i], nnop[i], 0]; + Tminus = Simplify @ Total @ Map[isospinminus[#, nnop2[#]]&, basisops]; + bzvc2bzop @ transformQStoIS[fnc @ basisops, Tminus] +]; + +quickISObasis[basisops__] := quickbasis[basisops, qsbasisvc]; +quickISOSZbasis[basisops__] := quickbasis[basisops, qszbasisvc]; +quickSU2basis[basisops__] := quickbasis[basisops, qbasisvc]; + +(* Product basis with double (isospin, etc.) symmetry *) +quickDBL[basisops1_, basisops2_, fnc_] := Module[{bz1, bz2}, + bz1 = fnc @ basisops1; + bz2 = fnc @ basisops2; + makebasis[Join[basisops1, basisops2]]; + (* Drop the second index *) + bz1 = mergebasis @ Map[{{#[[1, 1]]}, #[[2]]} &, bz1]; + bz2 = mergebasis @ Map[{{#[[1, 1]]}, #[[2]]} &, bz2]; + directproductbasis[bz1, bz2, Join[#1, #2] &] +]; + +quickDBLSZ[basisops1_, basisops2_, fnc_] := Module[{bz1, bz2}, + bz1 = fnc @ basisops1; + bz2 = fnc @ basisops2; + makebasis[Join[basisops1, basisops2]]; + mergebasis @ directproductbasis[bz1, bz2, {#1[[1]], #2[[1]], #1[[2]]+#2[[2]]} &] +]; + +(* bzQ[b] returns True if the argument b is of the correct form to + represent a basis: a list of {{quantum numbers}, {states}} pairs. *) +bzQ[{ {{__},{___}}.. }] = True; +bzQ[_] = False; + +(* makematricesbzop[] and makematricesbzvc[] produce a table of operator + matrices for all invariant subspaces. The functions are appropriate for + operators that are diagonal with respect to the conserved quantum numbers + that are used to enumerate the invariant subspaces, i.e. singlets with + respect to the group of symmetry operations which are explicitly taken + into account in the problem. For general operators, use + makeallmatricesbzop[] and makeallmatricesbzvc[] which combine all pairs + of invariant subspaces. *) + +makematricesbzop[op_, basis_?bzQ] := Map[ + {First[#], Simplify @ matrixrepresentationop[op, #[[2]], #[[2]]]}&, basis ]; + +makematricesbzvc[op_, basis_?bzQ] := Map[ + {First[#], Simplify @ matrixrepresentationvc[op, #[[2]], #[[2]]]}&, basis ]; + +(* Helper function for makeallmatricesbzvc : process a pair of invariant + subspaces. FNC is either matrixrepresentationvc or + matrixrepresentationop. *) + +mambpair[{qn1_, vecs1_}, {qn2_, vecs2_}, op_, FNC_] := Module[{pair, matrep}, + pair = {qn1, qn2}; + matrep = FNC[op, vecs1, vecs2]; + (* TRICK: Drop combinations of subspaces that do not contribute. *) + tot = Total[Abs[matrep], 2]; + If[tot =!= 0, {pair, matrep}, HoldComplete[Sequence[]]] +]; + +listofpairs[l1_, l2_] := Flatten[Outer[List, l1, l2, 1], 1]; + +(* A generalization of makematricesbzvc[] and makematricesbzop[] which + handle operators which are not diagonal in the quantum numbers that + enumerate the invariant subspaces. *) + +makeallmatricesbzvc[op_, basis_] := ReleaseHold @ + mambpair[Sequence @@ #, op, matrixrepresentationvc]& /@ + listofpairs[basis, basis]; + +makeallmatricesbzop[op_, basis_] := ReleaseHold @ + mambpair[Sequence @@ #, op, matrixrepresentationop]& /@ + listofpairs[basis, basis]; + +(* Used in maskconstants[] and other functions. *) +getallconstants[] := Union[listrealconstants, listcomplexconstants, + listintegerconstants, listfreeindexes]; + +(* maskconstants[] builds a list of replacement rules for all constants that +have been defined using snegrealconstants[] and similar functions. It +returns the rules and the inverse rules. This is useful, since +Eigensystem[] has problems when parameters in the matrix are defined to have +special properties. *) + +maskconstants[] := Module[{l}, + l = getallconstants[]; + l = Map[# -> ToExpression[ToString[#] <> ToString[Unique[]]] &, l]; + {l, Reverse /@ l} +]; + +(* Perform operation op while all known constants have been masked. *) + +maskOp[op_] := Module[{map, invmap, tmp, res}, + {map, invmap} = maskconstants[]; + tmp = Hold[op] /. map; + res = ReleaseHold[tmp]; + res /. invmap +]; + +(* Sometimes arguments needs to be evaluated before being passed to +a function. *) + +maskOp[op_, args__] := Module[{map, invmap, tmp, res}, + {map, invmap} = maskconstants[]; + tmp = {args} /. map; + res = ReleaseHold[op] @@ tmp; + res /. invmap +]; + +SetAttributes[maskOp, HoldFirst]; + +(* Cross product *) + +mcross[a_List, b_List] := { + nc[a[[2]], b[[3]]] - nc[a[[3]], b[[2]]], + nc[a[[3]], b[[1]]] - nc[a[[1]], b[[3]]], + nc[a[[1]], b[[2]]] - nc[a[[2]], b[[1]]] +}; + + +(* +++ Simplification rules +++ *) + +allsnegsimplifyrules = {Automatic, SnegSimplifyNumber, SnegSimplifySpin, + SnegSimplifyHop, SnegSimplifyHubbard, SnegSimplifyMisc}; + +SnegSimplifyNumber[expr_] := + expr //. {nc[op_[CR, j___, sigma_], op_[AN, j___, sigma_]] :> + HoldForm[number[op[j], sigma]] /; operatorQ[op], + Expand[z_. (HoldForm[number[i_, 1]] + HoldForm[number[i_, 0]]) ] :> + z HoldForm[number[i]]}; + +ruleSnegSimplifySpin = + Block[{op1, op2, rule1, rule2, rule3, l1, l2, ir1, ir2, i1, i2}, + (* Note: By setting upvalues, the following rules will be removed + when op1 and op2 go out of scope. *) + operatorQ[op1] ^= True; + operatorQ[op2] ^= True; + fermionQ[op1] ^= True; + fermionQ[op2] ^= True; + l1 = {HoldForm @ spinx[op1[i1]], + HoldForm @ spiny[op1[i1]], + HoldForm @ spinz[op1[i1]]}; + l2 = {HoldForm @ spinx[op2[i2]], + HoldForm @ spiny[op2[i2]], + HoldForm @ spinz[op2[i2]]}; + ir1 = {op1 -> a_?operatorQ, op2 -> b_?operatorQ}; + ir2 = {op1 -> a, op2 -> b}; + + (* Rules for a single spin operator *) + + rule1 = Thread[(spinxyz[op1[i1___]] /. ir1) -> (l1 /. ir2)]; + rule1a = Thread[(Expand[spinxyz[op1[i1___]]] /. ir1) -> (l1 /. ir2)]; + rule1b = Thread[(Expand[2spinxyz[op1[i1___]]] /. ir1) -> (2l1 /. ir2)]; + rule1c = Thread[(Expand[I spinxyz[op1[i1___]]] /. ir1) -> (I l1 /. ir2)]; + rule1d = Thread[(Expand[2I spinxyz[op1[i1___]]] /. ir1) -> (2I l1 /. ir2)]; + rule1e = Thread[(Simplify[2I spinxyz[op1[i1___]]] /. ir1) -> (2I l1 /. ir2)]; + rule1f = Thread[(Simplify[-2I spinxyz[op1[i1___]]] /. ir1) -> (-2I l1 /. ir2)]; + + (* Rules for products of spin operators *) + + rule2 = Thread[ + Flatten[Expand[ + 4 outer[spinxyz[op1[i1___]], spinxyz[op2[i2___]]]] /. + ir1, 1] -> + + Flatten[Expand[4 outer[l1, l2]] /. ir2, 1]]; + + rule2b = Thread[ + Flatten[Expand[ + outer[spinxyz[op1[i1___]], spinxyz[op2[i2___]]]] /. + ir1, 1] -> + + Flatten[Expand[outer[l1, l2]] /. ir2, 1]]; + + (* Some symmetrized combinations of spin - spin products *) + rule2a = { + (* xx + yy *) + Rule @@ (rule2[[1]] + rule2[[5]] /. + HoldPattern[i_ -> j_] :> Expand[-1/2{z_. i, z j}]), + (* xy + yx *) + Rule @@ (rule2[[2]] + rule2[[4]] /. + HoldPattern[i_ -> j_] :> Expand[I/2{z_. i, z j}]), + (* xz + zx *) + Rule @@ (rule2[[3]] + rule2[[7]] /. + HoldPattern[i_ -> j_] :> Expand[{z_. i, z j}]), + (* yz + zy *) + Rule @@ (rule2[[6]] + rule2[[8]] /. + HoldPattern[i_ -> j_] :> Expand[-I{z_. i, z j}]) + }; + + (* Rules for combining the components of spin operator pairs *) + rule3 = {Expand[z_. inner[l1, l2] /. + Join[ir1, {i1 -> i1___, i2 -> i2___}]] -> + (z HoldForm[spinspin[op1[i1], op2[i2]]] /. ir2)}; + + Join[rule3, + rule1, rule1a, rule1b, rule1c, rule1d, rule1e, rule1f, + rule2, rule2b, rule2a] + ]; + +SnegSimplifySpin[expr_] := Simplify[expr, TransformationFunctions -> + {Automatic, Function[x, x //. ruleSnegSimplifySpin]}]; + +ruleSnegSimplifyHop = Block[{op1, op2, ir1, ir2, rule1, rule2}, + operatorQ[op1] ^= True; + operatorQ[op2] ^= True; + fermionQ[op1] ^= True; + fermionQ[op2] ^= True; + ir1 = {op1 -> a_?operatorQ, op2 -> b_?operatorQ}; + ir2 = {op1 -> a, op2 -> b}; + rule1 = Expand[z_. hop[op1[i___], op2[j___]] /. ir1] -> + (z HoldForm[hop[op1[i], op2[j]]] /. ir2); + rule2 = Expand[z_. twohop[op1[i___], op2[j___]] /. ir1] -> + (z HoldForm[twohop[op1[i], op2[j]]] /. ir2); + {rule1, rule2} + ]; + +SnegSimplifyHop[expr_] := expr //. ruleSnegSimplifyHop; + +ruleSnegSimplifyHubbard = Block[{op1}, + operatorQ[op1] ^= True; + fermionQ[op1] ^= True; + (-hubbard[op1[i___]] /. {op1 -> a_?operatorQ}) -> -HoldForm[ + hubbard[a[i]]] + ]; + +SnegSimplifyHubbard[expr_] := expr //. ruleSnegSimplifyHubbard; + +ruleSnegSimplifyMisc = Block[{op1, op2, ir1, ir2, rule1, rule2}, + operatorQ[op1] ^= True; + operatorQ[op2] ^= True; + fermionQ[op1] ^= True; + fermionQ[op2] ^= True; + ir1 = {op1 -> a_?operatorQ, op2 -> b_?operatorQ}; + ir2 = {op1 -> a, op2 -> b}; + (* Sz.Sz + n.n rule *) + rule1 = (Expand[z_. (2spinz[op1[i___]]~nc~spinz[op2[j___]] - + 1/2chargecharge[op1[i___], op2[j___]]) ] /. ir1) -> + z (2nc[HoldForm[spinz[op1[i]]], HoldForm[spinz[op2[j]]]] - + 1/2HoldForm[chargecharge[op1[i], op2[j]]] /. ir2); + rule2 = (Expand[z_. (-2spinz[op1[i___]]~nc~spinz[op2[j___]] - + 1/2chargecharge[op1[i___], op2[j___]]) ] /. ir1) -> + z (-2nc[HoldForm[spinz[op1[i]]], HoldForm[spinz[op2[j]]]] - + 1/2HoldForm[chargecharge[op1[i], op2[j]]] /. ir2); + {rule1, rule2} +]; + +SnegSimplifyMisc[expr_] := expr //. ruleSnegSimplifyMisc; + +(* SnegSimplify[] works well for "simple" expressions *) + +SnegSimplify[expr_] := Module[{x}, + x = Collect[expr, getallconstants[] ]; + x = SnegSimplifyNumber[x]; + x = SnegSimplifyHubbard[x]; + x = SnegSimplifyHop[x]; + x = SnegSimplifySpin[x]; + x = SnegSimplifyMisc[x]; + x = Collect[x, getallconstants[] ] +]; + +SnegFullSimplify[expr_] := Simplify[expr, + TransformationFunctions -> { Automatic, SnegSimplifyNumber, + SnegSimplifyHubbard, SnegSimplifyHop, SnegSimplifySpin, + SnegSimplifyMisc, SnegSimplify }]; + +(* snegSeries[f, expr, n] computes the power series of function f +around 0 to order n. The result is returned in terms of powers +of expr. *) +snegSeries[f_, expr_, n_:5] := Module[{s, x}, + s = Normal @ Series[f[x], {x, 0, n}]; + s /. x^pw_. -> fastpow[expr, pw] (* pow -> fastpow, 25.8.2010 *) +]; + +snegSeries[f_[expr_], n_:5] := snegSeries[f, expr, n]; +snegSeries[E^expr_, n_:5] := snegSeries[Exp, expr, n]; + +(* Example: snegSeries[Exp[c]] *) + +(* SNTR[in_] generates a test line to be used in test - sneg.m for + regression testing. *) + +SetAttributes[SNTR, HoldFirst]; +SNTR[in_] := Module[{out}, + out = Release[in]; + fn[expr_] := ToString @ InputForm @ expr; + strip[str_String] := StringDrop[StringDrop[str, -1], 9]; + "test[" <> strip @ fn @HoldForm[in] <> ", " <> fn @out <> "];"]; + +(* snegAssuming evaluates expr for assumed value of an argument. + TagSet must be used to define the assumption. + Upvalues and downvalues of the argument are saved and restored + after the evaluation. *) +(* snegAssumingHoldForm is similar to snegAssuming, but it returns + the result enwrapped by HoldForm[]. *) + +SetAttributes[snegAssuming, HoldAllComplete]; +SetAttributes[snegAssumingHoldForm, HoldAllComplete]; + + +snegAssumingHoldForm[x:HoldPattern[{TagSet[_,__]..}], expr_] := Module[ +{storeup, storedown, res, vars}, + vars = ReleaseHold[Hold[x] /. HoldPattern[TagSet[var_, defs__]] :> var]; + storeup = Map[{#, UpValues[#]}&, vars]; + storedown = Map[{#, DownValues[#]}&, vars]; + ReleaseHold[x]; + res = ReleaseHold[expr]; + res = Apply[HoldForm, {res}]; + Map[(UpValues[Evaluate @ First[#]] = Last[#])&, storeup]; + Map[(DownValues[Evaluate @ First[#]] = Last[#])&, storedown]; + res +]; + +snegAssumingHoldForm[HoldPattern[x:TagSet[var_, def__]], expr_] := + snegAssumingHoldForm[{x}, expr]; + +snegAssuming[x___] := snegAssumingHoldForm[x] /. HoldForm[a_] -> a; + +(* Simplify KroneckerDelta and UnitStep expressions. *) +SimplifyKDfunc[expr_] := (expr //. { + KroneckerDelta[k_, l_] UnitStep[k_] UnitStep[-l_] -> 0, + KroneckerDelta[k_, l_] UnitStep[-k_] UnitStep[l_] -> 0, + UnitStep[k_] UnitStep[-k_] -> 0, + UnitStep[k_]+UnitStep[-k_] -> 1, + 1 - UnitStep[-k_] -> UnitStep[k], + UnitStep[-k_] -1 -> -UnitStep[k], + HoldPattern[KroneckerDelta[x__]^n_Integer] -> KroneckerDelta[x], + HoldPattern[UnitStep[x__]^n_Integer] -> UnitStep[x] + }); + +SimplifyKD[expr_] := Simplify[expr, +TransformationFunctions -> {Automatic, SimplifyKDfunc}]; + +FullSimplifyKD[expr_] := FullSimplify[expr, +TransformationFunctions -> {Automatic, SimplifyKDfunc}]; + +(* "Thread" spin projection values enlisted in l over all operators +that appear in expr. *) +snegSpinThread[expr_, l_List] := + Plus @@ Map[expr /. o_?operatorQ[i___] :> o[i, #] &, l]; + +sum2list[HoldPattern[Plus[i___]]] := {i}; +sum2list[i_] := {i}; + +(* Similar to matrixrepresentationvc[], but this function +returns a sparse matrix and it is significantly faster. *) +(* UNTESTED! *) +matrixrepresentationvcsparse[a_, l_List] := Module[{}, + SparseArray[Select[Flatten[ + Table[ + (sum2list @ Collect[ap[a, l[[i]]], l]) /. + x_. v_vc :> ({i, Position[l, v][[1, 1]]} -> x), + {i, Length[l]}], + 1], (#=!=0)&]] +]; + +(* Similar to matrixrepresentationvc[], but significantly faster. *) +(* NOTE: only applicable for simple monomial vectors!! *) +matrixrepresentationvcfast[a_, l_List] := Module[{n = Length[l]}, + Map[ + indexvalue2list[ + Sort[ + (sum2list @ Collect[ap[a, #], l]) /. + x_. v_vc :> {Position[l, v][[1,1]], x} + ], + n] &, + l] // Transpose +]; + +matrixrepresentationvcfast[a_, l1_List, l2_List] := Module[{n = Length[l1]}, + Map[ + indexvalue2list[ + Sort[ + (sum2list @ Collect[ap[a, #], l1]) /. + x_. v_vc :> Module[{pos = Position[l1, v]}, + If[pos==={},0,{pos[[1,1]], x}] ] + ], + n] &, + l2] // Transpose +]; + +(* Transform a list of {index, value} pairs into a vector with given +values at indexed positions. *) +indexvalue2list[iv_List, n_Integer] := + Fold[ReplacePart[#1, Last[#2], First[#2]] &, Table[0, {n}], + Select[iv, (# =!= 0)&] ]; + +(* Added 25.3.2011 *) +indexvalue2list[{0}, n_Integer]:= Table[0, {n}]; +indexvalue2list[{0.}, n_Integer]:= Table[0, {n}]; +indexvalue2list[{0.+0. I}, n_Integer]:= Table[0, {n}]; + +(* Symmetrized spin-spin correlation for a list of operators. *) +spinspinsymmetric[ops_List] := Module[{l}, + l=spinspin @@@ Subsets[ops, {2}]; + Total[l]/Length[l] // Simplify +]; + +(* Show Hamiltonian matrices in an easy-to-comprehend formatted form. *) +showbzmat[bz_?bzQ] := + MatrixForm[Map[{First[#], MatrixForm @ Last[#]}&, bz]]; + +(* Make a dispatch table which allows very fast computation of operator.vc operations. *) +(* vc basis must be previously declared by makebasis[], so that BASIS contains the +relevant operators which appear in the operator. *) +opdispatch[op_] := Module[{states, resultstates}, + states = vc @@@ Tuples[{0, 1}, Length[BASIS]]; + resultstates = nc[op, states]; + Dispatch[Thread[states -> resultstates]] +]; + +(* Routines for producing tensor product basis for SU(2) multiplets *) + +qscombinefnc[s_, {q1_, _}, {q2_, _}] := {q1 + q2, s}; +qsspincombinefnc[s_, {q1_, _}, {___}] := {q1, s}; + +SU2basistensorproduct[l1_, l2_, + spinfnc1_, spinfnc2_, + combineqnfnc_, + fixspin1fnc_, fixspin2fnc_] := + mergebasis @ + Flatten[Outer[ + SU2merge[#1, #2, spinfnc1, spinfnc2, + combineqnfnc, fixspin1fnc, fixspin2fnc] &, + l1, l2, 1], + 2]; + +SU2merge[{qn1_, l1_}, {qn2_, l2_}, + spinfnc1_, spinfnc2_, + combineqnfn_, + fixspin1fnc_, fixspin2fnc_] := + Module[{s1, s2, states}, + s1 = spinfnc1[qn1]; + s2 = spinfnc2[qn2]; + states = + Partition[ + Flatten[ + Outer[ + SU2combine[#1, #2, {s1, s2}, fixspin1fnc, fixspin2fnc] &, + l1, l2] + ], + 2]; + states = mergebasis @ + Map[{combineqnfn[First[#], qn1, qn2], {Last[#]}} &, states] + ]; + +SU2combine[op1_, op2_, {s1_, s2_}, fixspin1fnc_, fixspin2fnc_] := + Table[ + {s, Sum[ + With[{cg = ClebschGordan[{s1, sz1}, {s2, sz2}, {s, s}]}, + If[cg =!= 0, + cg nc[fixspin1fnc[op1, s1, sz1] , fixspin2fnc[op2, s2, sz2]], + 0] + ], {sz1, -s1, s1}, {sz2, -s2, s2}]}, + {s, Abs[s1 - s2], s1 + s2} + ]; + +fixspin[op_, s_, s_] := op; +fixspin[op_, s_, sz_] /; sz < s := Nest[spindown, op, s - sz]; + +fixspinket[k_, s_, s_] := k; +fixspinket[k_, s_, sz_] /; sz < s := k /. ket[i_] :> ket[i - (s - sz)]; + +(* Miscellaneous useful routines *) + +(* Returns a n-by-n matrix of zeros *) +ZeroMatrix[n_] := 0 IdentityMatrix[n]; + +(* SU(3) symmetry *) +(* The defining representation: Gell-Mann matrices. *) +su3g[1] = {{0, 1, 0}, {1, 0, 0}, {0, 0, 0}}; +su3g[2] = {{0, -I, 0}, {I, 0, 0}, {0, 0, 0}}; +su3g[3] = {{1, 0, 0}, {0, -1, 0}, {0, 0, 0}}; +su3g[4] = {{0, 0, 1}, {0, 0, 0}, {1, 0, 0}}; +su3g[5] = {{0, 0, -I}, {0, 0, 0}, {I, 0, 0}}; +su3g[6] = {{0, 0, 0}, {0, 0, 1}, {0, 1, 0}}; +su3g[7] = {{0, 0, 0}, {0, 0, -I}, {0, I, 0}}; +su3g[8] = 1/Sqrt[3] {{1, 0, 0}, {0, 1, 0}, {0, 0, -2}}; +su3t[i_] = su3g[i]/2; + +(* Construct all combination from a list of alternatives *) + +allcombinations[l_List] := allcombinations[l, {{}}]; +allcombinations[{}, g_List] := g; +allcombinations[l_List, g_List] := + allcombinations[Rest[l], Flatten[Outer[Append, g, First[l], 1], 1]]; + + +(* New much improved implementation of the VEV calculation using +Wick's theorem. *) + +sneglinearoperatorFirst[vevwicknew]; + +vevwicknew[ nc[___, x:op_?fermionQ[__]] ] /; isannihilation[x] := 0; +vevwicknew[ nc[x:op_?fermionQ[__], ___] ] /; iscreation[x] := 0; + +vevwicknew[op_?operatorQ[__]] := 0; (* single operator *) +vevwicknew[HoldPattern[l_nc]] /; OddQ[Length[l]] := 0; + +vevwicknew[HoldPattern[l_nc]] /; EvenQ[Length[l]] := + Module[{oplist, l1, l2, crops, conjcrops, pos, l3, contributions}, + oplist = List @@ l; + l1 = indecesCR[l]; + l2 = indecesAN[l]; + + (* Extract creation operators and conjugate them *) + crops = Map[l[[#]] &, l1]; + conjcrops = conj /@ crops; + + (* Locate these in the operator list *) + pos = Map[Position[oplist, #, {1}] &, conjcrops]; + pos = Replace[pos, {n_} :> n, {2}]; + + (* Generate all allowed permutations *) + l3 = allcombinations[pos]; + l3 = Select[l3, DuplicateFreeQ]; + + (* Calculate the contractions *) + contributions = Map[contract[oplist, Transpose[{l1, #}]] &, l3]; + + Total[contributions /. {dd[_]->0}] +]; + +End[]; (* End Private Context *) + +EndPackage[]; + diff --git a/TBpack/TBpack.m b/TBpack/TBpack.m index a650feb..78216c1 100644 --- a/TBpack/TBpack.m +++ b/TBpack/TBpack.m @@ -6,13 +6,13 @@ (* :Title: TBpack *) (* :Author: Vasil A. Saroka <40.ovasil@gmail.com> *) (* :Context: TBpack` *) -(* :Version: 0.2.1 *) -(* :Date: 2020-03-18 *) +(* :Version: 0.3.0 *) +(* :Date: 2021-03-04 *) (* :Mathematica Version: 10.0+ *) (* :Copyright: (c) 2020 Vasil A. Saroka *) -BeginPackage["TBpack`", {"TBpack`UnitcellGenerators`", "TBpack`DataAnalysis`","TBpack`Optics`"}]; +BeginPackage["TBpack`", {"TBpack`UnitcellGenerators`", "TBpack`DataAnalysis`","TBpack`Optics`","TBpack`Sneg`"}]; Unprotect[Evaluate[$Context<>"*"]]; (* taken from CustomTicks package *)