# # !!!!!! The package is currently under revision - not all commands are here # # Andriy Zhugayevych azh@ukr.net # LatticeTools package (this under revision) # created 25.05.2004 modified - see version below # ################################################################################ #root: LatticeTools #hfl: LatticeTools #toc: _Overview LatticeTools:=module() option package; export ModuleLoad, Setup, gfunRd, propRd, # Miscellaneous functions BuildW, BuildA, BuildH, # Builders GreensF, StatSol, ZeroExpan, # Solvers PathExpansion, # Analyzers Velocity, DiffusionTensor, DiffusionLength; # Infinite periodic lattices local version,t,s,x,v; ModuleLoad:=proc() version:=20230928; # BasicTools package is required if not(member('BasicTools',packages())) then try with(BasicTools) catch: error "Cannot find the required BasicTools package: %1",lastexception end end; # Default settings Setup(); NULL end: #hfl: Setup #toc: _Setup #LatticeTools[Setup] Setup:=proc({printout::boolean:=false}) ProcessSetupArgs([_rest],[],['version']) end: ################################################################################ #cat: Miscellaneous functions #hfl: Functions gfunRd:=(d,s,r)->(2*Pi)^(-d/2)*(sqrt(s)/r)^(d/2-1)*BesselK(d/2-1,sqrt(s)*r): #hfl: Functions propRd:=(d,t,r)->(4*Pi*t)^(-d/2)*exp(-r^2/4/t): ################################################################################ #cat: Builders #hfl: BuildW BuildW:=proc(n::posint,s::{"Matrix","RandomMatrix","Generate","Sample","GamEps"},gen,{sym::boolean:=false}) local W,Gam,Eps,i,j; if (s="Matrix") then W:=Matrix(n,gen,`if`(sym,shape=symmetric,NULL),_rest) elif (s="RandomMatrix") then W:=LinearAlgebra[RandomMatrix](n,'generator'=gen,`if`(sym,shape=symmetric,NULL),_rest) elif (s="Generate") then W:=Matrix(n,RandomTools[Generate](gen,'makeproc'=true),`if`(sym,shape=symmetric,NULL),_rest) elif (s="Sample") then W:=Matrix(n,convert(Statistics[Sample](gen,n^2),list),_rest); if sym then W:=Matrix(W,shape=symmetric,_rest) end elif (s="GamEps") then Gam:=Matrix(n,convert(Statistics[Sample](gen[1],n^2),list)); Eps:=Statistics[Sample](gen[2],n); W:=Matrix(n,(i,j)->`if`(i=j,0,exp(Eps[i])*Gam[i,j]),_rest) end; for i from 1 to n do W[i,i]:=0 end; for i from 1 to n do for j from 1 to n do if (W[i,j]<0) then error "W[%1,%2]=%3<0",i,j,W[i,j] end end end; if (s="GamEps") then W,Gam,Eps else W end end: #hfl: BuildA BuildA:=proc(W::Matrix,U::Vector:=Vector(op(1,W)[1],0),$) local n,i,j,k; n:=op(1,W)[1]; Matrix(n,(i,j)->`if`(i=j,-U[i]-add(W[i,k],k=1..n),W[i,j]),shape=MatrixOptions(W,shape),datatype=MatrixOptions(W,datatype)) end: #hfl: BuildH BuildH:=proc(tbH0::table,k::{name,indexable}:=undefined,{symmetric::boolean:=false,printout::boolean:=false},$) local tbH,mmts,d,n,mmt,mmt2,rng,o,N,mts,m,mt,n2,mtsid,m2,H,Hk,vars; tbH:=copy(tbH0); mmts:=[indices(tbH)]; d:=nops(mmts[1])-2; n:=max(seq(op(..2,mmt),mmt=mmts)); if symmetric then for mmt in mmts do mmt2:=[mmt[2],mmt[1],op(-mmt[3..])]; if member(mmt2,mmts) then if (tbH[op(mmt2)]<>tbH[op(mmt)]) then error("Nonsymmetric tbH%1=%2<>%3=tbH%4",mmt2,tbH[op(mmt2)],tbH[op(mmt)],mmt) end else tbH[op(mmt2)]:=tbH[op(mmt)] end end; mmts:=[indices(tbH)] end; rng:=[seq([(min,max)(map2(op,2+o,mmts))],o=1..d)]; N:=1+max(map(abs,map(op,rng))); mts:=Sort(convert({seq(mmt[2..],mmt=mmts)} minus {seq([m,0$d],m=1..n)},list),mt->((add(mt[1+o]^2,o=1..d)*N^d+add(N^(o-1)*abs(mt[1+o]),o=1..d))*2^d+add(2^(o-1)*`if`(mt[1+o]<0,1,0),o=1..d))*n+mt[1]); mts:=[seq([m,0$d],m=1..n),op(mts)]; n2:=nops(mts); if printout then printf("d=%d, n=%d, n2=%d, |tbH|=%d, rng=%s\n",d,n,n2,nops(mmts),StringTools[DeleteSpace](sprintf("%a",rng))) end; mtsid:=table([seq(mts[m2]=m2,m2=1..n2)]); H,Hk:=Matrix(n,n2,`if`(hastype([entries(tbH,'nolist')],float),datatype=float,NULL)),Matrix(n); for mmt in mmts do H[mmt[1],mtsid[mmt[2..]]]:=tbH[op(mmt)]; Hk[op(..2,mmt)]:=Hk[op(..2,mmt)]+tbH[op(mmt)]*exp(I*add(k[o]*mmt[2+o],o=1..d)) end; for m from 1 to n do Hk[m,m]:=evalc(Hk[m,m]) end; vars:=indets(H); if (vars<>{}) then Hk:=map(collect,Hk,vars); if printout then printf("vars=%s\n",StringTools[DeleteSpace](sprintf("%a",vars))) end end; H,Hk,mts,op(mtsid) end: ################################################################################ #cat: Solvers #hfl: GreensF GreensF:=proc(s,W::Matrix,U::Vector:=Vector(op(1,W)[1],0),$) local n,Gi; n:=op(1,W)[1]; Gi:=Matrix(n,(i,j)->`if`(i=j,s+U[i]+add(W[i,k],k=1..n),-W[i,j]),shape=MatrixOptions(W,shape),datatype=`if`(type(s,float),float,MatrixOptions(W,datatype))); Gi^(-1) end: #hfl: StatSol StatSol:=proc(W::Matrix,i00::nonnegint:=0,simplifyf::procedure:=simplify,{printout::boolean:=false},$) local A,params,i0,pi1,C,pi,i,v; A:=BuildA(W); params:=select(type,indets(A),name); i0:=`if`(i00=0,`if`(params=[],op(MinIdx([seq(A[i,i],i=1..op(1,A)[1])])),1),i00); pi1:=LinearAlgebra[LinearSolve](LinearAlgebra[Transpose](A[[..i0-1,i0+1..]$2]),LinearAlgebra[Transpose](-A[i0,[..i0-1,i0+1..]])); C:=`minus`(select(type,indets(pi1),name),params); if (C<>{}) then pi1:=simplify(subs(seq(v=0,v=C),pi1)) end; pi:=Vector[row](op(1,pi1)+1,`if`(VectorOptions(pi1,'datatype')=float[8],'datatype'=float,NULL)); pi[i0]:=1; pi[[..i0-1,i0+1..]]:=pi1; pi:=LinearAlgebra[Normalize](pi,1); if (params<>{} or simplifyf<>simplify) then pi:=simplifyf(pi) end; if (printout and params={}) then printf("Norm(pi.A)= %.1e (L1), %.1e (Linf)\n",LinearAlgebra[Norm](pi.A,1),LinearAlgebra[Norm](pi.A,infinity)) end; pi end: #hfl: StatSol ZeroExpan:=proc(W::Matrix,i0::nonnegint:=0,simplifyf::procedure:=simplify,{printout::boolean:=false},$) local dtdecl,pi,n,T,Ti,i,lsi,W2,pi2,R2,R,j,v; dtdecl:=`if`(MatrixOptions(W,'datatype')=float[8],'datatype'=float,NULL); pi:=StatSol(W,i0,simplifyf,':-printout'=printout); n:=op(1,pi); T:=Matrix(n,n-1,(i,j)->`if`(i=j+1,1,0)-pi[j+1],dtdecl); Ti:=Matrix(n-1,n,dtdecl); for i from 1 to n-1 do Ti[i,1]:=-1; Ti[i,i+1]:=1 end; try pi,T.LinearAlgebra[MatrixInverse](-Ti.BuildA(W).T).Ti catch: lsi:=[seq(`if`(pi[i]=0,NULL,i),i=1..n)]; if printout then printf("Singular matrix, will use subspace [%{c,}d%s]\n",Vector(lsi[..min(16,nops(lsi))]),`if`(nops(lsi)>16,",...","")) end; W2:=W[lsi,lsi]; if (add(abs(v),v=W2)=0) then return pi,Matrix(n,dtdecl) end; pi2,R2:=ZeroExpan(W2,':-printout'=printout); R:=Matrix(n,dtdecl); R[lsi,lsi]:=R2; pi,R end end: ################################################################################ #cat: Analyzers #hfl: PathExpansion PathExpansion:=proc(H::Matrix,E::numeric,oi::posint,of::posint,vmin::numeric,Nmax::posint,subsys::list(posint):=[oi,of], {Lmax::posint:=99,neglect::numeric:=1e-3,printout::boolean:=false,digits::posint:=3,nprint::posint:=9},$) local n,env,vmin2,maxorder2,tau,T,Lmax2,L,vmax,ls,o,o1,o2,v; n:=Dim2(H)[1]; env:=remove(member,[$1..n],subsys); vmin2:=neglect*vmin; maxorder2:=-round(log10(vmin))-digits; if printout then printf("oi=%d, of=%d, Ei=%.*f, Ef=%.*f, E=%.*f\n",oi,of,digits,H[oi,oi],digits,H[of,of],digits,E) end; tau:=Matrix(n,(o1,o2)->`if`(o1=o2,0,H[o1,o2]/(E-H[o2,o2])),datatype=float); T:=table([0=[[1,oi]]]); Lmax2:=Lmax; for L from 1 to Lmax2 do T[L]:=remove(v->abs(v[1])abs(v[1]),T[L])),digits,'maxorder'=maxorder2)) end; if (nops(T[L])>Nmax) then T[L]:=Sort(T[L],v->-abs(v[1]))[..Nmax]; WARNING("Too many terms, cut at %1",sprintf("%.0e",abs(T[L][-1][1]))) end end; if printout then printf("Capping and removing small terms:\n") end; for L from 0 to Lmax2 do T[L]:=remove(v->abs(v[1])[v[1]*H[v[-1],of],op(3..,v)],T[L])) end; for L from Lmax2 by -1 to 1 while (T[L]=[]) do unassign('T[L]') end; Lmax2:=L; if printout then for L from 0 to Lmax2 do printf("%2d %8d %*.*f %s\n",L,nops(T[L]),digits+4,digits,add(v[1],v=T[L]),FormatFloat(`if`(T[L]=[],0,max(map(v->abs(v[1]),T[L]))),digits,'maxorder'=maxorder2)) end; printf("Hbare=%.*f, Hrenorm=%.*f, Hpath=%.*f\nLargest elements:\n",digits,H[oi,of],digits,H[oi,of]+H[oi,env].(E-H[env,env])^(-1).H[env,of],digits,add(add(v[1],v=T[L]),L=0..Lmax2)); vmax:=max(seq(seq(abs(v[1]),v=T[L]),L=0..Lmax2)); ls:=Sort([seq(op(remove(v->abs(v[1])-abs(v[1])); seq(printf("%*.*f [%{c,}d]\n",digits+4,digits,v[1],Vector(v[2..])),v=ls[..min(nprint,nops(ls))]) end; op(T) end: ################################################################################ #cat: Infinite periodic lattices #hfl: Velocity Velocity:=proc( sites::Array, W::Matrix, M::Matrix:=LinearAlgebra[IdentityMatrix](Dim2(sites)[2]), replicas::{'undefined',Vector}:='undefined', i0::nonnegint:=0, simplifyf::procedure:=simplify, { printout::boolean:=false },$) local n,n2,d,dtdecl,repl,i,pi; n,n2:=Dim2(W); d:=Dim2(sites)[2]; dtdecl:=`if`(MatrixOptions(W,datatype)=float or MatrixOptions(W,datatype)[0]=float,datatype=float,NULL); for i from 1 to n do W[i,i]:=0 end: if (replicas='undefined') then repl:=Vector(n,i->[]): for i from 1 to n2 do repl[sites[i,0]]:=[op(repl[sites[i,0]]),i] end else repl:=replicas end; pi:=StatSol(Matrix(n,(i,j)->`if`(i=j,0,add(W[i,l],l=repl[j])),dtdecl),i0,simplifyf,':-printout'=printout); Vector[column](pi.Matrix(n,d,(i,o)->add(W[i,j]*sites[j,o],j=1..n2),dtdecl).LinearAlgebra[Transpose](M),dtdecl) end: #hfl: Velocity DiffusionTensor:=proc( sites::Array, W::Matrix, M::Matrix:=LinearAlgebra[IdentityMatrix](Dim2(sites)[2]), replicas::{'undefined',Vector}:='undefined', i0::nonnegint:=0, simplifyf::procedure:=simplify, { printout::boolean:=false },$) local n,n2,d,dtdecl,repl,i,Wp,Wpp,pi,R,uv,j,l,o,o1,o2; n,n2:=Dim2(W); d:=Dim2(sites)[2]; dtdecl:=`if`(MatrixOptions(W,datatype)=float or MatrixOptions(W,datatype)[0]=float,datatype=float,NULL); for i from 1 to n do W[i,i]:=0 end: if (replicas='undefined') then repl:=Vector(n,i->[]): for i from 1 to n2 do repl[sites[i,0]]:=[op(repl[sites[i,0]]),i] end else repl:=replicas end; pi,R:=ZeroExpan(Matrix(n,(i,j)->`if`(i=j,0,add(W[i,l],l=repl[j])),dtdecl),i0,simplifyf,':-printout'=printout); Wp:=Vector(d,o->Matrix(n,(i,j)->add(W[i,l]*sites[l,o],l=repl[j]),dtdecl)); Wpp:=Matrix(d,(o1,o2)->Matrix(n,(i,j)->add(W[i,l]*sites[l,o1]*sites[l,o2],l=repl[j]),dtdecl)); uv:=Vector(n,1): Matrix(M.Matrix(d,(o1,o2)->(pi.Wpp[o1,o2].uv+pi.Wp[o1].R.Wp[o2].uv+pi.Wp[o2].R.Wp[o1].uv)/2,shape=symmetric,dtdecl).LinearAlgebra[Transpose](M),shape=symmetric,dtdecl) end: #hfl: Velocity DiffusionLength:=proc( sites::Array, W::Matrix, U::Vector, points::Matrix, M::Matrix:=LinearAlgebra[IdentityMatrix](Dim2(sites)[2]), replicas::{'undefined',Vector}:='undefined', { printout::boolean:=false },$) local n,n2,d,dtdecl,repl,i,G,Wp,Wpp,Up,Upp,uv; n,n2:=Dim2(W); d:=Dim2(sites)[2]; dtdecl:=`if`(MatrixOptions(W,datatype)=float or MatrixOptions(W,datatype)[0]=float,datatype=float,NULL); for i from 1 to n do W[i,i]:=0 end: if (replicas='undefined') then repl:=Vector(n,i->[]): for i from 1 to n2 do repl[sites[i,0]]:=[op(repl[sites[i,0]]),i] end else repl:=replicas end; G:=GreensF(0,Matrix(n,(i,j)->`if`(i=j,0,add(W[i,l],l=repl[j])),dtdecl),U); Wp:=Vector(d,o->Matrix(n,(i,j)->add(W[i,l]*sites[l,o],l=repl[j]),dtdecl)); Wpp:=Matrix(d,(o1,o2)->Matrix(n,(i,j)->add(W[i,l]*sites[l,o1]*sites[l,o2],l=repl[j]),dtdecl)); Up:=Vector(d,o->Vector(n,i->U[i]*points[o,i],dtdecl)); Upp:=Matrix(d,(o1,o2)->Vector(n,i->U[i]*points[o1,i]*points[o2,i],dtdecl)); uv:=Vector(n,1): Vector(n,i->Matrix(M.Matrix(d,(o1,o2)->G[i,..].(Wpp[o1,o2].uv+Wp[o1].G.Wp[o2].uv+Wp[o2].G.Wp[o1].uv+Wp[o1].G.Up[o2]+Wp[o2].G.Up[o1]+Upp[o1,o2])/2,shape=symmetric,dtdecl).LinearAlgebra[Transpose](M),shape=symmetric,dtdecl)) end: ModuleLoad() end module: