###################################################### # # Fock space Hamiltonian truncation for QED in 1+1 dimensions # Author: Ian Lim (itlim@ucdavis.edu), adapted from work by Ryzhkov and Vitale # July 2020 # ###################################################### import scipy from scipy import pi import numpy as np import scipy.sparse.linalg import scipy.sparse import scipy.interpolate from operator import attrgetter from math import factorial from phi1234 import Matrix from statefuncs import Basis, NotInBasis, omega, State from oscillators import NormalOrderedOperator as NOO from qedops import uspinor, vspinor, FermionOperator from qedstatefuncs import FermionBasis, FermionState from dressedstatefuncs import DressedFermionBasis, DressedFermionState from dressedstatefuncs import ZeroModeRaisingOperator, ZeroModeLoweringOperator from dressedstatefuncs import calculateAxialCharge import collections import renorm import itertools from progress.bar import ChargingBar import time tol = 0.0001 """ P denotes spatial parity, while K field parity. For now only the P-even sector is implemented """ class Schwinger(): """ main class Attributes ---------- L : float Circumference of the circle on which to quantize m : float Mass of the scalar field Emax: float Maximum energy cutoff All basic dictionaries in this class have two keys, 1 and -1, corresponding to values of k-parity. h0: dict of Matrix objects Values are Matrix objects corresponding to the free Hamiltonian potential: dict of dict of Matrix Values are dictionaries where the keys are orders in the field (0,2,4 correspond to phi^0,phi^2,phi^4) and the values are the potential matrices at each order stored as Matrix objects When actually diagonalizing the Hamiltonian, we further restrict to a subspace with some different nmax, possibly? H: dict of Matrix objects Like h0, but on a restricted subspace defined by basis rather than fullBasis V: dict of dict of Matrix Like potential but restricted to the basis subspace """ def __init__(self): self.L = None self.m = None self.Emax = None self.h0 = None self.potential = None self.h0OffDiag = None self.h0Sub = None self.H = None self.V = None self.h0OffDiagSub = None self.eigenvalues = None self.eigsrenlocal = None self.eigsrensubl = None self.eigenvectors = None self.basis = None self.fullBasis = None def buildFullBasis(self,L,m,Emax,bcs="periodic"): """ Builds the full Hilbert space basis """ self.L=float(L) self.m=float(m) #self.fullBasis = FermionBasis(self.L, Emax, self.m, bcs=bcs) self.fullBasis = DressedFermionBasis(self.L, Emax, self.m, bcs=bcs) # TO-DO: test this method. We only have one interaction term so it's just V def buildBasis(self,Emax): """ Builds the Hilbert space basis for which the Hamiltonian to actually diagonalize is calculated (in general it's a subspace of fullBasis) Note that this is called in phi4eigs, but not in generating the original potential matrix and free Hamiltonian. """ self.basis = DressedFermionBasis(self.L, Emax, self.m, nmax=self.fullBasis.nmax, bcs=self.fullBasis.bcs) # We use the vector length (nmax) of the full basis. # In this way we can compare elements between the two bases self.Emax = float(Emax) self.V = self.potential.sub(self.basis, self.basis).M.tocoo() self.h0Sub = self.h0.sub(self.basis,self.basis).M.tocoo() self.h0OffDiagSub = self.h0OffDiag.sub(self.basis,self.basis).M.tocoo() def buildMatrix(self): """ Builds the full Hamiltonian in the basis of the free Hamiltonian eigenvectors. This is computationally intensive. It can be skipped by loading the matrix from file """ """ Possible speedups: Simplify the diagonal operator loops? Can we apply a numpy mask to simplify the loops without the conditionals? Basically this loops over the operators at each order in phi and stores all the normal order operators explicitly in lists. """ L=self.L m=self.m basis = self.fullBasis lookupBasis = self.fullBasis Emax = basis.Emax nmax = basis.nmax """ We split operators into diagonal and off-diagonal operators. The idea is that for off-diagonal operators, we can just compute the image of half the operators and generate the conjugates by taking the transpose. That is, if an operator X takes state A to B, then X^\dagger should take B to A. """ diagOps = {0: None, 2:None, 4:None} offdiagOps = {0: None, 2:None, 4:None} """ diagOps[0] = [ NOO([],[],L,m) ] offdiagOps[0] = [] diagOps[2] = [ NOO([a],[a],L,m, extracoeff=2.) for a in range(-nmax,nmax+1) ] offdiagOps[2] = [ NOO([a,-a],[],L,m,extracoeff=comb(a,-a)) for a in range(-nmax,nmax+1) if a<=-a<=nmax and omega(a,L,m)+omega(-a,L,m) <= Emax+tol] diagOps[4] = [ NOO([a,b],[c,a+b-c],L,m, extracoeff=6.*comb(a,b)*comb(c,a+b-c)) for a in range(-nmax,nmax+1) for b in range (a,nmax+1) for c in range(-nmax,nmax+1) if ( c<=a+b-c<=nmax and (a,b) == (c,a+b-c) and -Emax-tol <= omega(a,L,m)+omega(b,L,m) - omega(c,L,m)-omega(a+b-c,L,m) <=Emax+tol)] offdiagOps[4] = [ NOO([a,b,c,-a-b-c],[],L,m,extracoeff=comb(a,b,c,-a-b-c)) for a in range(-nmax,nmax+1) for b in range (a,nmax+1) for c in range(b,nmax+1) if c<=-a-b-c<=nmax and omega(a,L,m)+omega(b,L,m) + omega(c,L,m)+omega(-a-b-c,L,m)<= Emax+tol] \ + [ NOO([a,b,c],[a+b+c],L,m, extracoeff = 4. * comb(a,b,c)) for a in range(-nmax, nmax+1) for b in range (a,nmax+1) for c in range(b,nmax+1) if (-nmax<=a+b+c<=nmax and -Emax-tol <= omega(a,L,m)+omega(b,L,m)+ omega(c,L,m)-omega(a+b+c,L,m) <=Emax+tol)] \ + [ NOO([a,b],[c,a+b-c],L,m, extracoeff = 6. * comb(a,b)*comb(c,a+b-c)) for a in range(-nmax,nmax+1) for b in range (a,nmax+1) for c in range(-nmax,nmax+1) if ( c<=a+b-c<=nmax and (a,b) != (c,a+b-c) and sorted([abs(a),abs(b)]) < sorted([abs(c),abs(a+b-c)]) and -Emax-tol <= omega(a,L,m)+omega(b,L,m)- omega(c,L,m)-omega(a+b-c,L,m) <=Emax+tol)] """ #save as h0 a lookupBasis.size by 1 sparse matrix initialized to zeros #really just a single column #and store the bases as h0[k].basisI and h0[k].basisJ #self.h0[k] = Matrix(lookupBasis, basis) tempEnergies = np.empty(basis.size) #tempEnergies = np.array([basis[j].energy # for j in range(basis.size)]) for j in range(basis.size): tempEnergies[j] = basis[j].energy #self.h0[k].finalize() temph0 = scipy.sparse.diags(tempEnergies,format="coo") self.h0 = Matrix(lookupBasis,basis,temph0) """ We build the potential in this basis. self.potential is a dictionary with two keys corresponding to k=1 and k=-1. Each of the entries is a list of sparse matrices. Each sparse matrix consists of a sum of off-diagonal components and diagonal components in COO format. Possible speedup: could we apply each operator to all states in one shot? Then we would just loop over operators. Right now, for each state, we apply each operator one at a time and mark down its image in the corresponding column (the origin state) and row (the image state). If we apply an operator to a vector of states, we will get a vector of their images which can be collectively checked for validity by checking the dictionary keys/energy ranges, and then this is just the matrix form of the operator. """ #psidaggerpsi = self.generateOperators() # can change this to enumerate if we want to specialize to diagonal # operators # the use of only the upper triangular matrix and then just transposing # is a great way to save time # maybe later should just enumerate all 16 ops by hand? Seems like # a pain though. allops = self.generateOperators3(nmax) potential = Matrix(lookupBasis,basis) bar = ChargingBar("States", max=basis.size) for state in basis: newcolumn = np.zeros(lookupBasis.size) for op in allops: try: normalization, index = op.apply2(basis,state) if index != None: newcolumn[index] += normalization except NotInBasis: pass ''' for k in np.arange(-nmax,nmax+1): for op2 in psidaggerpsi[-k]: n, newstate= op2._transformState(state,returnCoeff=True) if n == 0: continue for op1 in psidaggerpsi[k]: try: normalization, index = op1.apply2(basis,newstate) if (index != None): #for ease of comparison we can put the 1/2L later newcolumn[index] += (normalization * n / (2*self.L*k**2)) except NotInBasis: pass ''' potential.addColumn(newcolumn) bar.next() bar.finish() potential.finalize() #print(potential.M.toarray()) isSymmetric = (np.array_equal(potential.M.toarray(), potential.M.toarray().T)) assert isSymmetric, "Matrix not symmetric (hermitian)" # if isSymmetric: # print("Matrix is symmetric") # # if np.any(np.diag(potential.M.toarray())): # print("nonzero diagonal entries") self.potential = potential #off diagonal h0 elements h0offdiag = Matrix(lookupBasis,basis) raising_op = ZeroModeRaisingOperator() lowering_op = ZeroModeLoweringOperator() for state in basis: newcolumn = np.zeros(lookupBasis.size) axialcharge = calculateAxialCharge(state) if axialcharge == 0: potential.addColumn(newcolumn) pass for op in [raising_op,lowering_op]: try: normalization, index = op.apply2(basis,state) if index != None: #this second factor is because h0 'counts' the number #of fermion insertions and adds or subtracts factors #of A1. possibly there's a minus sign here? newcolumn[index] += normalization * (-1*axialcharge) except NotInBasis: pass h0offdiag.addColumn(newcolumn) h0offdiag.finalize() #print(potential.M.toarray()) isSymmetric = (np.array_equal(h0offdiag.M.toarray(), h0offdiag.M.toarray().T)) assert isSymmetric, "Matrix not symmetric (hermitian)" # if isSymmetric: # print("Matrix is symmetric") # # if np.any(np.diag(potential.M.toarray())): # print("nonzero diagonal entries") #include normalization factors of 1/sqrt(2)*1/pi^1/4 self.h0OffDiag = h0offdiag *(2**-0.5)*(pi**-0.25) # for each order (0,2,4) in phi """ for n in offdiagOps.keys(): offdiag_V = Matrix(lookupBasis, basis) diagonal = np.zeros(basis.size) # for each state in the basis for j in range(basis.size): newcolumn = np.zeros(lookupBasis.size) # for each off-diagonal operator at a given order for op in offdiagOps[n]: try: # apply this operator to find whether the # new state is still in the basis (x,i) = op.apply(basis,j,lookupBasis) # if so, add the corresponding value to the matrix # this is basically writing the effects of the # operator in the basis of the free states if(i != None): newcolumn[i]+=x except NotInBasis: pass offdiag_V.addColumn(newcolumn) # for each diagonal operator at the same order n for op in diagOps[n]: (x,i) = op.apply(basis,j,lookupBasis) # It should be j=i if i!= None: if i != j: raise RuntimeError('Non-diagonal operator') diagonal[i]+=x offdiag_V.finalize() diag_V = scipy.sparse.spdiags(diagonal,0,basis.size,basis.size) self.potential[n] = (offdiag_V+offdiag_V.transpose()+Matrix(lookupBasis, basis, diag_V)).to('coo')*self.L """ # TO-DO: update for QED def saveMatrix(self, fname): """ Saves the free Hamiltonian and potential matrices to file """ t = (fname, self.L, self.m, \ self.fullBasis.Emax, self.fullBasis.nmax, self.fullBasis.bcs, \ self.h0.M.data,self.h0.M.row,self.h0.M.col, \ self.potential.M.data,self.potential.M.row,self.potential.M.col, \ ) scipy.savez(*t) def loadMatrix(self, fname): """ Loads the free Hamiltonian and potential matrices from file """ f = scipy.load(fname) self.L = f['arr_0'].item() self.m = f['arr_1'].item() Emax = f['arr_2'].item() nmax = f['arr_3'].item() bcs = f['arr_4'].item() self.buildFullBasis(L=self.L, m=self.m, Emax=Emax, bcs=bcs) basisI = self.fullBasis basisJ = self.fullBasis print(basisI.size) self.h0 = Matrix(basisI, basisJ, scipy.sparse.coo_matrix((f['arr_5'], (f['arr_6'], f['arr_7'])), shape=(basisI.size, basisJ.size))) self.potential = Matrix(basisI, basisJ, scipy.sparse.coo_matrix((f['arr_8'], (f['arr_9'], f['arr_10'])), shape=(basisI.size, basisJ.size))) def generateOperators(self): """ Generates a dictionary of normal-ordered operators corresponding to the operator (\psi^\dagger \psi)_k for all values of k within range. Returns ------- opsList : dict of FermionOperators A dictionary whose keys are values of k, i.e. which Fourier mode of psi^\dagger \psi is under consideration, and whose entries are the corresponding set of operators. """ opsList = {} # we will generate (psi^\dagger psi)_k for k in np.arange(-self.fullBasis.nmax,self.fullBasis.nmax+1): opsList[k] = [] if k == 0: continue # note: what is the range on n? revisit. All valid ns such that # -nmax <= k+n <= nmax and -nmax <= n <= nmax for n in np.arange(-self.fullBasis.nmax,self.fullBasis.nmax+1): if not (-self.fullBasis.nmax <= k-n <= self.fullBasis.nmax): continue # zero mode spinor wavefunctions fixed, can remove this? # if k-n == 0 or n == 0: # continue coeff = np.vdot(uspinor(n-k,self.L,self.m,normed=True), uspinor(n,self.L,self.m,normed=True)) adaggera = FermionOperator([n-k],[n],[],[],self.L,self.m, extracoeff=coeff/np.sqrt(self.L), normed=True) #n.b. -1 from anticommutation of bdagger and adagger coeff = -1 * np.vdot(uspinor(n-k,self.L,self.m,normed=True), vspinor(-n,self.L,self.m,normed=True)) bdaggeradagger = FermionOperator([n-k],[],[-n],[],self.L,self.m, extracoeff=coeff/np.sqrt(self.L), normed=True) coeff = np.vdot(vspinor(k-n,self.L,self.m,normed=True), uspinor(n,self.L,self.m,normed=True)) ba = FermionOperator([],[n],[],[k-n],self.L,self.m, extracoeff=coeff/np.sqrt(self.L),normed=True) # -1 from anticommuting b and b dagger coeff = -1 * np.vdot(vspinor(k-n,self.L,self.m,normed=True), vspinor(-n,self.L,self.m,normed=True)) bdaggerb = FermionOperator([],[],[-n],[k-n],self.L,self.m, extracoeff=coeff/np.sqrt(self.L),normed=True) #the anticommutator is always trivial because k-n = -n -> k=0 #and we've precluded this possiblity since k != 0 opsList[k] += [adaggera, bdaggeradagger, ba, bdaggerb] return opsList def udotu(self,k1,k2): return np.vdot(uspinor(k1,self.L,self.m,normed=True), uspinor(k2,self.L,self.m,normed=True)) def udotv(self,k1,k2): return np.vdot(uspinor(k1,self.L,self.m,normed=True), vspinor(k2,self.L,self.m,normed=True)) def vdotu(self,k1,k2): return np.vdot(vspinor(k1,self.L,self.m,normed=True), uspinor(k2,self.L,self.m,normed=True)) def vdotv(self,k1,k2): return np.vdot(vspinor(k1,self.L,self.m,normed=True), vspinor(k2,self.L,self.m,normed=True)) def makeInteractionOps(self,clist,dlist,anticlist,antidlist,nmax, spinors,coeff=1.,deltacondition=None): """ Takes a set of field labels for (anti)particle creation/annihilation operators and generates all corresponding strings of operators with momenta up to nmax, subject to (optional) delta function constraints. Parameters ---------- clist : list of int labels for which fields the particle creation operators come from, e.g. [1,3] indicates that the first and third fields contributed a daggers. dlist : list of int labels for which fields the particle annihilation operators come from anticlist : list of int labels for which fields the antiparticle creation operators come from antidlist : list of int labels for which fields the antiparticle annihilation operators come from nmax : int maximum value of wavenumber in the truncation spinors : string A four-character string representing the contracted spinor wavefunctions to be computed. For example, udagger u vdagger v is just given as "uuvv". coeff : float, optional Extra coefficients for this operator. The default is 1. deltacondition : tuple, optional Any pairs of momenta that are to be set equal by Kronecker deltas, e.g. (1,2) is \delta_{k1,k2}. The default is None. Returns ------- ops : list of FermionOperators All fermion operators matching the input parameters. The normalization is for 1/2k^2 * |\psi^\dagger \psi|^2, i.e. it contains the spinor inner products and also the 1/L factor. The g^2 then multiplies the matrix entries of V. """ momenta = np.array([[k1,k2,k3,-k1-k2-k3] for k1 in np.arange(-nmax,nmax+1) for k2 in np.arange(-nmax,nmax+1) for k3 in np.arange(-nmax,nmax+1) if k1 + k2 != 0 and abs(k1+k2+k3) <= nmax ]) # note: easily modified for multiple delta functions # just make it a list of tuples and iterate over the masking conditions # a minus sign is needed for the delta conditions with this convention # since we have a_k1 and a^\dagger_{-k2} if deltacondition: mask = momenta[:,deltacondition[0]-1] == -momenta[:,deltacondition[1]-1] momenta = momenta[mask] # for kvals in momenta: # assert(kvals[2]+kvals[3] != 0) # assert(kvals[0]+kvals[1] + kvals[2] + kvals[3] == 0) # return the right functions for the spinor inner products # with this convention, vs and udaggers get minus signs spinor_lookup = {"uu" : lambda k1,k2: self.udotu(-k1,k2), "uv" : lambda k1,k2: self.udotv(-k1,-k2), "vu" : lambda k1,k2: self.vdotu(k1,k2), "vv" : lambda k1,k2: self.vdotv(k1,-k2) } assert len(spinors) == 4 firstspinor = spinor_lookup[spinors[0:2]] secondspinor = spinor_lookup[spinors[2:4]] ops = [] for kvals in momenta: spinorfactor = (firstspinor(kvals[0],kvals[1]) * secondspinor(kvals[2],kvals[3])) if spinorfactor == 0: continue # print(firstspinor(kvals[0],kvals[1])) # print(secondspinor(kvals[2],kvals[3])) # print(spinorfactor) # print(kvals) ksquared = 1/(kvals[0]+kvals[1])**2 #creation operators get the minus sign on k ops += [FermionOperator(-kvals[np.array(clist,dtype=int)-1], kvals[np.array(dlist,dtype=int)-1], -kvals[np.array(anticlist,dtype=int)-1], kvals[np.array(antidlist,dtype=int)-1], self.L,self.m, extracoeff=coeff*spinorfactor*ksquared/(2*self.L), normed=True)] return ops def generateOperators2(self,nmax): """ Parameters ---------- nmax : TYPE DESCRIPTION. Returns ------- None. """ """ operators are named by their non-normal ordered versions and how many basic operators are in the final version, e.g. this first operator is the normal ordered version of adagger a adagger a, which looks like adagger adagger a a and there is also a 2-operator term adagger a from the one anticommutation """ ops = [] adagger_a_adagger_a_4 = self.makeInteractionOps([1,3],[2,4],[],[], nmax=nmax, coeff=-1., spinors="uuuu") ops += adagger_a_adagger_a_4 adagger_a_adagger_a_2 = self.makeInteractionOps([1],[4],[],[], nmax=nmax, deltacondition=(2,3), spinors="uuuu") ops += adagger_a_adagger_a_2 adagger_a_adagger_bdagger_4 = self.makeInteractionOps([1,3],[2],[4],[], nmax=nmax, spinors="uuuv") ops += adagger_a_adagger_bdagger_4 adagger_a_adagger_bdagger_2 = self.makeInteractionOps([1],[],[4],[], nmax=nmax, coeff=-1., deltacondition=(2,3), spinors="uuuv") ops += adagger_a_adagger_bdagger_2 adagger_bdagger_adaggger_a_4 = self.makeInteractionOps([1,3],[4],[2],[], nmax=nmax, coeff=-1., spinors="uvuu") ops += adagger_bdagger_adaggger_a_4 adagger_a_b_bdagger_4 = self.makeInteractionOps([1],[2],[4],[3], nmax=nmax, spinors="uuvv") ops += adagger_a_b_bdagger_4 ''' adagger_a_b_bdagger_2 = self.makeInteractionOps([1],[2],[],[], nmax=nmax, deltacondition=(3,4), spinors="uuvv") ops += adagger_a_b_bdagger_2 ''' adagger_bdagger_b_bdagger_4 = self.makeInteractionOps([1],[],[2,4],[3], nmax=nmax, coeff=-1., spinors="uvvv") ops += adagger_bdagger_b_bdagger_4 ''' adagger_bdagger_b_bdagger_2 = self.makeInteractionOps([1],[],[2],[], nmax=nmax, coeff=-1., deltacondition=(3,4), spinors="uvvv") ops += adagger_bdagger_b_bdagger_2 ''' adagger_bdagger_adagger_bdagger_4 = self.makeInteractionOps([1,3],[],[2,4],[], nmax=nmax, coeff=-1., spinors="uvuv") ops += adagger_bdagger_adagger_bdagger_4 b_bdagger_adagger_bdagger_4 = self.makeInteractionOps([3],[],[2,4],[1], nmax=nmax, spinors="vvuv") ops += b_bdagger_adagger_bdagger_4 b_bdagger_adagger_bdagger_2_1 = self.makeInteractionOps([3],[],[2],[], nmax=nmax, deltacondition=(1,4), spinors="vvuv") ops += b_bdagger_adagger_bdagger_2_1 ''' b_bdagger_adagger_bdagger_2_2 = self.makeInteractionOps([3],[],[4],[], nmax=nmax, coeff=-1., deltacondition=(1,2), spinors="vvuv") ops += b_bdagger_adagger_bdagger_2_2 ''' adagger_bdagger_b_a_4 = self.makeInteractionOps([1],[4],[2],[3], nmax=nmax, coeff=-1., spinors="uvvu") ops += adagger_bdagger_b_a_4 b_a_adagger_bdagger_4 = self.makeInteractionOps([3],[2],[4],[1], nmax=nmax, coeff=-1., spinors="vuuv") ops += b_a_adagger_bdagger_4 b_a_adagger_bdagger_2_1 = self.makeInteractionOps([3],[2],[],[], nmax=nmax, coeff=-1., deltacondition=(1,4), spinors="vuuv") ops += b_a_adagger_bdagger_2_1 b_a_adagger_bdagger_2_2 = self.makeInteractionOps([],[],[4],[1], nmax=nmax, coeff=-1., deltacondition=(2,3), spinors="vuuv") ops += b_a_adagger_bdagger_2_2 # note: there is a total delta function term here but it just shifts # the vacuum energy so we omit it. It would be delta_14 delta_23 vuuv. b_bdagger_b_bdagger_4 = self.makeInteractionOps([],[],[2,4],[1,3], nmax=nmax, coeff=-1., spinors="vvvv") ops += b_bdagger_b_bdagger_4 b_bdagger_b_bdagger_2_1 = self.makeInteractionOps([],[],[2],[3], nmax=nmax, deltacondition=(1,4), spinors="vvvv") ops += b_bdagger_b_bdagger_2_1 ''' b_bdagger_b_bdagger_2_2 = self.makeInteractionOps([],[],[2],[1], nmax=nmax, coeff=-1., deltacondition=(3,4), spinors="vvvv") ops += b_bdagger_b_bdagger_2_2 b_bdagger_b_bdagger_2_3 = self.makeInteractionOps([],[],[4],[3], nmax=nmax, coeff=-1., deltacondition=(1,2), spinors="vvvv") ops += b_bdagger_b_bdagger_2_3 ''' #next the hermitian conjugate terms #we can probably do this with the off diagonal trick later b_a_adagger_a_4 = self.makeInteractionOps([3],[2,4],[],[1], nmax=nmax, spinors="vuuu") ops += b_a_adagger_a_4 b_a_adagger_a_2 = self.makeInteractionOps([],[4],[],[1], nmax=nmax, deltacondition=(2,3), spinors="vuuu") ops += b_a_adagger_a_2 adagger_a_b_a_4 = self.makeInteractionOps([1],[2,4],[],[3], nmax=nmax, coeff=-1., spinors="uuvu") ops += adagger_a_b_a_4 b_bdagger_adagger_a_4 = self.makeInteractionOps([3],[4],[2],[1], nmax=nmax, spinors="vvuu") ops += b_bdagger_adagger_a_4 ''' b_bdagger_adagger_a_2 = self.makeInteractionOps([3],[4],[],[], nmax=nmax, deltacondition=(1,2), spinors="vvuu") ops += b_bdagger_adagger_a_2 ''' b_bdagger_b_a_4 = self.makeInteractionOps([],[4],[2],[1,3], nmax=nmax, coeff=-1., spinors="vvvu") ops += b_bdagger_b_a_4 ''' b_bdagger_b_a_2 = self.makeInteractionOps([],[4],[],[3], nmax=nmax, deltacondition=(1,2), spinors="vvvu") ops += b_bdagger_b_a_2 ''' b_a_b_a_4 = self.makeInteractionOps([],[2,4],[],[1,3], nmax=nmax, coeff=-1., spinors="vuvu") ops += b_a_b_a_4 b_a_b_bdagger_4 = self.makeInteractionOps([],[2],[4],[1,3], nmax=nmax, spinors="vuvv") ops += b_a_b_bdagger_4 ''' b_a_b_bdagger_2_1 = self.makeInteractionOps([],[2],[],[1], nmax=nmax, deltacondition=(3,4), spinors="vuvv") ops += b_a_b_bdagger_2_1 ''' b_a_b_bdagger_2_2 = self.makeInteractionOps([],[2],[],[3], nmax=nmax, coeff=-1., deltacondition=(1,4), spinors="vuvv") ops += b_a_b_bdagger_2_2 return ops def generateOperators3(self,nmax): """ Parameters ---------- nmax : TYPE DESCRIPTION. Returns ------- None. """ """ operators are named by their non-normal ordered versions and how many basic operators are in the final version, e.g. this first operator is the normal ordered version of adagger a adagger a, which looks like adagger adagger a a and there is also a 2-operator term adagger a from the one anticommutation """ ops = [] adagger_a_adagger_a_4 = self.makeInteractionOps([1,3],[2,4],[],[], nmax=nmax, coeff=-1., spinors="uuuu") ops += adagger_a_adagger_a_4 adagger_a_adagger_a_2 = self.makeInteractionOps([1],[4],[],[], nmax=nmax, deltacondition=(2,3), spinors="uuuu") ops += adagger_a_adagger_a_2 adagger_a_adagger_bdagger_4 = self.makeInteractionOps([1,3],[2],[4],[], nmax=nmax, coeff=2., spinors="uuuv") ops += adagger_a_adagger_bdagger_4 adagger_a_adagger_bdagger_2 = self.makeInteractionOps([1],[],[4],[], nmax=nmax, coeff=-1., deltacondition=(2,3), spinors="uuuv") ops += adagger_a_adagger_bdagger_2 # adagger_bdagger_adaggger_a_4 = self.makeInteractionOps([1,3],[4],[2],[], # nmax=nmax, # coeff=-1., # spinors="uvuu") # ops += adagger_bdagger_adaggger_a_4 adagger_a_b_bdagger_4 = self.makeInteractionOps([1],[2],[4],[3], nmax=nmax, coeff=2., spinors="uuvv") ops += adagger_a_b_bdagger_4 adagger_bdagger_b_bdagger_4 = self.makeInteractionOps([1],[],[2,4],[3], nmax=nmax, coeff=-2., spinors="uvvv") ops += adagger_bdagger_b_bdagger_4 adagger_bdagger_adagger_bdagger_4 = self.makeInteractionOps([1,3],[],[2,4],[], nmax=nmax, coeff=-1., spinors="uvuv") ops += adagger_bdagger_adagger_bdagger_4 # b_bdagger_adagger_bdagger_4 = self.makeInteractionOps([3],[],[2,4],[1], # nmax=nmax, # spinors="vvuv") # ops += b_bdagger_adagger_bdagger_4 b_bdagger_adagger_bdagger_2_1 = self.makeInteractionOps([3],[],[2],[], nmax=nmax, deltacondition=(1,4), spinors="vvuv") ops += b_bdagger_adagger_bdagger_2_1 adagger_bdagger_b_a_4 = self.makeInteractionOps([1],[4],[2],[3], nmax=nmax, coeff=-2., spinors="uvvu") ops += adagger_bdagger_b_a_4 # b_a_adagger_bdagger_4 = self.makeInteractionOps([3],[2],[4],[1], # nmax=nmax, # coeff=-1., # spinors="vuuv") # ops += b_a_adagger_bdagger_4 b_a_adagger_bdagger_2_1 = self.makeInteractionOps([3],[2],[],[], nmax=nmax, coeff=-1., deltacondition=(1,4), spinors="vuuv") ops += b_a_adagger_bdagger_2_1 b_a_adagger_bdagger_2_2 = self.makeInteractionOps([],[],[4],[1], nmax=nmax, coeff=-1., deltacondition=(2,3), spinors="vuuv") ops += b_a_adagger_bdagger_2_2 # note: there is a total delta function term here but it just shifts # the vacuum energy so we omit it. It would be delta_14 delta_23 vuuv. b_bdagger_b_bdagger_4 = self.makeInteractionOps([],[],[2,4],[1,3], nmax=nmax, coeff=-1., spinors="vvvv") ops += b_bdagger_b_bdagger_4 b_bdagger_b_bdagger_2_1 = self.makeInteractionOps([],[],[2],[3], nmax=nmax, deltacondition=(1,4), spinors="vvvv") ops += b_bdagger_b_bdagger_2_1 #next the hermitian conjugate terms #we can probably do this with the off diagonal trick later b_a_adagger_a_4 = self.makeInteractionOps([3],[2,4],[],[1], nmax=nmax, coeff=2., spinors="vuuu") ops += b_a_adagger_a_4 b_a_adagger_a_2 = self.makeInteractionOps([],[4],[],[1], nmax=nmax, deltacondition=(2,3), spinors="vuuu") ops += b_a_adagger_a_2 # adagger_a_b_a_4 = self.makeInteractionOps([1],[2,4],[],[3], # nmax=nmax, # coeff=-1., # spinors="uuvu") # ops += adagger_a_b_a_4 # b_bdagger_adagger_a_4 = self.makeInteractionOps([3],[4],[2],[1], # nmax=nmax, # spinors="vvuu") # ops += b_bdagger_adagger_a_4 b_bdagger_b_a_4 = self.makeInteractionOps([],[4],[2],[1,3], nmax=nmax, coeff=-2., spinors="vvvu") ops += b_bdagger_b_a_4 b_a_b_a_4 = self.makeInteractionOps([],[2,4],[],[1,3], nmax=nmax, coeff=-1., spinors="vuvu") ops += b_a_b_a_4 # b_a_b_bdagger_4 = self.makeInteractionOps([],[2],[4],[1,3], # nmax=nmax, # spinors="vuvv") # ops += b_a_b_bdagger_4 b_a_b_bdagger_2_2 = self.makeInteractionOps([],[2],[],[3], nmax=nmax, coeff=-1., deltacondition=(1,4), spinors="vuvv") ops += b_a_b_bdagger_2_2 return ops def setcouplings(self, g): self.g = float(g) def renlocal(self,Er): self.g0r, self.g2r, self.g4r = renorm.renlocal(self.g2,self.g4,self.Emax,Er) self.Er = Er def computeHamiltonian(self, ren=False): """ Computes the (renormalized) Hamiltonian to diagonalize ren : if True, computes the eigenvalue with the "local" renormalization procedure, otherwise the "raw" eigenvalues """ # note: V (if using generateOperators2) is |psi^\dagger psi|^2/2k^2. # The 1/L factor is accounted for. self.H = self.h0Sub - self.V*self.g**2 + self.h0OffDiagSub*self.g**0.5 """ if not(ren): self.H[k] = self.h0Sub[k] + self.V[k][2]*self.g2 + self.V[k][4]*self.g4 else: self.H[k] = self.h0Sub[k] + self.V[k][0]*self.g0r + self.V[k][2]*self.g2r + self.V[k][4]*self.g4r """ def computeEigval(self, ren=False, corr=False, sigma=0, n=10): """ Diagonalizes the Hamiltonian and possibly computes the subleading renormalization corrections ren (bool): it should have the same value as the one passed to computeHamiltonian() corr (bool): if True, computes the subleading renormalization corrections, otherwise not. n (int): number of lowest eigenvalues to compute sigma : value around which the Lanczos method looks for the lowest eigenvalue. """ (self.eigenvalues, eigenvectorstranspose) = scipy.sparse.linalg.eigsh( self.H, k=n, sigma=sigma, which='LM', return_eigenvectors=True) """ if not ren: (self.eigenvalues[k], eigenvectorstranspose) = scipy.sparse.linalg.eigsh(self.H[k], k=n, sigma=sigma, which='LM', return_eigenvectors=True) else: (self.eigsrenlocal[k], eigenvectorstranspose) = scipy.sparse.linalg.eigsh(self.H[k], k=n, sigma=sigma, which='LM', return_eigenvectors=True) """ eigenvectors = eigenvectorstranspose.T # TO-DO: update this for QED if corr: print("Adding subleading corrections to k="+str(k), " eigenvalues") self.eigsrensubl[k] = np.zeros(n) cutoff = 5. for i in range(n): cbar = eigenvectors[i] if abs(sum([x*x for x in cbar])-1.0) > 10**(-13): raise RuntimeError('Eigenvector not normalized') Ebar = self.eigsrenlocal[k][i] self.eigsrensubl[k][i] += Ebar ktab, rentab = renorm.rensubl(self.g2, self.g4, Ebar, self.Emax, self.Er, cutoff=cutoff) tckren = { } tckren[0] = scipy.interpolate.interp1d(ktab,rentab.T[0],kind='linear') tckren[2] = scipy.interpolate.interp1d(ktab,rentab.T[1],kind='linear') tckren[4] = scipy.interpolate.interp1d(ktab,rentab.T[2],kind='linear') for nn in (0,2,4): #for a,b,Vab in itertools.izip(self.V[k][nn].row,self.V[k][nn].col,self.V[k][nn].data): for a,b,Vab in zip(self.V[k][nn].row,self.V[k][nn].col,self.V[k][nn].data): if a > b: continue elif a == b: c = 1 else: c = 2 Eab2= (self.basis[k][a].energy + self.basis[k][b].energy)/2. coeff = tckren[nn](Eab2) self.eigsrensubl[k][i] += c * coeff * cbar[a] * cbar[b] * Vab def vacuumE(self, ren="raw"): return self.eigenvalues[0] # implement others with corrections later """ if ren=="raw": return self.eigenvalues[1][0] elif ren=="renlocal": return self.eigsrenlocal[1][0] elif ren=="rensubl": return self.eigsrensubl[1][0] else: raise ValueError("Wrong argument") # The vacuum is K-even """ def spectrum(self, ren="raw"): if ren=="raw": eigs = self.eigenvalues elif ren=="renlocal": eigs = self.eigsrenlocal elif ren=="rensubl": eigs = self.eigsrensubl else: raise ValueError("Wrong argument") return scipy.array([x-self.vacuumE(ren=ren) for x in eigs[1:]]) """ # Now subtract vacuum energies if k==1: return scipy.array([x-self.vacuumE(ren=ren) for x in eigs[k][1:]]) elif k==-1: return scipy.array([x-self.vacuumE(ren=ren) for x in eigs[k]]) else: raise ValueError("Wrong argument") """