[KeldyshGF-svncheckins] [svn commit /var/www/svn/KeldyshGF] r268 - in trunk/kgf: . exact_diagonalization

svn at fysik.dtu.dk svn at fysik.dtu.dk
Fri Nov 16 15:10:50 CET 2012


Author: strange
Date: 2012-11-16 15:10:50 +0100 (Fri, 16 Nov 2012)
New Revision: 268

Added:
   trunk/kgf/exact_diagonalization/
   trunk/kgf/exact_diagonalization/benzene_ppp.py
   trunk/kgf/exact_diagonalization/exactd.py
   trunk/kgf/exact_diagonalization/exactd2.py
   trunk/kgf/exact_diagonalization/exactd3.py
   trunk/kgf/exact_diagonalization/exactdsparse.py
   trunk/kgf/exact_diagonalization/sparsetools.py
   trunk/kgf/exact_diagonalization/test.py
   trunk/kgf/exact_diagonalization/test_hubbard.py
   trunk/kgf/exact_diagonalization/test_hubbard_sparse.py
   trunk/kgf/exact_diagonalization/test_ppp_sparse.py
   trunk/kgf/exact_diagonalization/tmp.dat
Log:
tools for exact-diagonalization for el-el and el-ph interactions

Added: trunk/kgf/exact_diagonalization/benzene_ppp.py
===================================================================
--- trunk/kgf/exact_diagonalization/benzene_ppp.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/benzene_ppp.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,53 @@
+import numpy as np
+from numpy import linalg as npla
+from scipy import sparse as sp
+from scipy.sparse import linalg as spla
+import pylab as pl
+from exactdsparse import ExactDiag
+
+t = -2.539
+U = 10.06
+f = 14.397
+pos = 1.4 * np.asarray([[np.cos(a), np.sin(a), 0.0] 
+                         for a in np.arange(6)*2*np.pi/6.])
+n = len(pos)
+d2 = np.zeros((n, n))
+for i in range(n):
+    for j in range(n):
+        d2[i,j] = sum((pos[i]-pos[j])**2)
+
+d = np.sqrt(d2)
+h1 = t * ((d<1.5)*(d>0.0)).astype(float)
+h1[0, 0] += 1.0e-5 #break symmetry
+V = f / np.sqrt((f/U)**2 + d2)
+h_ion = np.zeros((n,n))
+for i in range(n):
+    h_ion[i,i] += -0.5 * V[i,i]
+    for j in range(n):
+        if j!=i: 
+            h_ion[i,i] += -1.0 * V[i,j]
+
+
+calc = ExactDiag(n)
+V2 = calc.get_twobody_PPP_ham(V)
+h2 = calc.get_onebody_ham(h1 + h_ion)
+H = h2 + V2
+print 'Lowest energy:', calc.get_lowest_energy(H)
+
+indices5 = calc.get_subspace_indices(npart=5, sz=0.5)
+indices6 = calc.get_subspace_indices(npart=6, sz=0.0)
+indices7 = calc.get_subspace_indices(npart=7, sz=0.5)
+
+calc.print_expectation_values(H, indices5)
+calc.print_expectation_values(H, indices6)
+calc.print_expectation_values(H, indices7)
+
+e7,sz7,s7,s27,n7 = calc.get_expectation_values(H, indices7)
+e6,sz6,s6,s26,n6 = calc.get_expectation_values(H, indices6)
+e5,sz5,s5,s25,n5 = calc.get_expectation_values(H, indices5)
+
+Egap = e5[0] + e7[0]-2 * e6[0]
+
+print 'Egap=%.3f' % Egap
+assert abs(Egap-11.3746) < 1.0e-3
+

Added: trunk/kgf/exact_diagonalization/exactd.py
===================================================================
--- trunk/kgf/exact_diagonalization/exactd.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/exactd.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,63 @@
+import numpy as np
+from numpy import ascontiguousarray as asc
+from numpy import linalg as la
+import pylab as pl
+
+def anticom(x,y):
+    return np.dot(x,y) + np.dot(y,x)
+
+# one single particle state
+P1 = np.identity(4,int)
+S1z = np.zeros((4,4),int)
+c1da = np.zeros((4,4),int)
+c1db = np.zeros((4,4),int)
+c1da[1,0] = 1
+c1da[3,2] = 1
+c1db[2,0] = 1
+c1db[3,1] = -1
+c1a = asc(c1da.T)
+c1b = asc(c1db.T)
+P1[1,1] = -1
+P1[2,2] = -1
+S1z[1,1] = 1
+S1z[2,2] = -1
+s1z = 0.5
+N1 = np.dot(c1da,c1a) + np.dot(c1db,c1b)
+
+# two single particle states
+c2da = np.kron(P1, c1da) # c_{2\uparrow}^\dagger
+c2db = np.kron(P1, c1db) # c_{2\downarrow}^\dagger
+
+c2a = asc(c2da.T)
+c2b = asc(c2db.T)
+
+c1da = np.kron(c1da, np.identity(4, int))
+c1db = np.kron(c1db, np.identity(4, int))
+c1a = asc(c1da.T)
+c1b = asc(c1db.T)
+
+N2 = np.kron(N1, np.identity(4,int)) + np.kron(np.identity(4,int), N1)
+P2 = np.kron(P1, P1) 
+S2z = np.kron(S1z, np.identity(4,int)) + np.kron(np.identity(4,int), S1z)
+
+print c1da
+print c1db
+print c2da
+print c2db
+# hamiltonian
+t = -2.0
+h = (np.dot(c1da, c2a) + np.dot(c1db, c2b) + np.dot(c2da, c1a) + 
+     np.dot(c2db, c1b)) * t
+
+def find_nparticle_subspace(npart, N):
+    return [i for i, n in enumerate(N.diagonal()) if n==npart]
+
+indices = find_nparticle_subspace(2, N2)
+ht = asc(h.take(indices,axis=0).take(indices,axis=1))
+#print np.around(ht,1)
+S2zt = asc(S2z.take(indices, axis=0).take(indices,axis=1))
+eigs, vecs = la.eig(ht)
+Sz_evaluated = np.around(np.dot(np.dot(vecs.T, S2zt*0.5), vecs),1).diagonal()
+#print N2.diagonal()
+#print Sz_evaluated
+#print np.around(eigs,1)

Added: trunk/kgf/exact_diagonalization/exactd2.py
===================================================================
--- trunk/kgf/exact_diagonalization/exactd2.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/exactd2.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,81 @@
+import numpy as np
+from numpy import ascontiguousarray as asc
+from numpy import linalg as la
+import pylab as pl
+
+# one single particle state
+P1 = np.identity(4,int)
+S1z = np.zeros((4,4),int)
+c1da = np.zeros((4,4),int)
+c1db = np.zeros((4,4),int)
+c1da[1,0] = 1
+c1da[3,2] = 1
+c1db[2,0] = 1
+c1db[3,1] = -1
+c1a = asc(c1da.T)
+c1b = asc(c1db.T)
+P1[1,1] = -1
+P1[2,2] = -1
+S1z[1,1] = 1
+S1z[2,2] = -1
+N1 = np.dot(c1da,c1a) + np.dot(c1db,c1b)
+
+def get_cds_and_cs(n):
+    cds = [[c1da, c1db]]
+    cs = [[c1a,c1b]]
+    Pop = P1.copy()
+    Nop = N1.copy()
+    for i in range(n-1):
+        I = np.identity(4**(i+1), int)
+        cda = np.kron(Pop, cds[i][0])
+        cdb = np.kron(Pop, cds[i][1])
+        ca = asc(cda.T)
+        cb = asc(cdb.T)
+
+        cds[i][0] = np.kron(cds[i][0], I)
+        cds[i][1] = np.kron(cds[i][1], I)
+        cs[i][0] = asc((cds[i][0]).T)
+        cs[i][1] = asc((cds[i][1]).T)
+        cds.append([cda,cdb])
+        cs.append([ca,cb])
+
+        Pop = np.kron(Pop, P1)
+        Nop = np.kron(Nop, np.identity(4,int)) + np.kron(I, N1)
+        return cds, cs, Nop
+
+def get_1body_ham(h, cds, cs):
+    n = h.shape[0]
+    H = np.zeros((4**n, 4**n), h.dtype) 
+    for i in range(n):
+        for j in range(n):
+            for s in range(2):
+                H += h[i,j] * np.dot(cds[i][s], cs[j][s])
+
+    return H
+
+def find_N_subspace(N, Nop):
+    return [i for i, n in enumerate(Nop.diagonal()) if n==N]
+
+def get_N_subspace(H, Nop, N):
+    indices = find_N_subspace(N, Nop)
+    return asc(H.take(indices, axis=0).take(indices, axis=1))
+
+# two basis functions
+#h = np.zeros((2,2), float)
+#h[0,1] = h[1,0] = -2.0
+#cds, cs, Nop = get_cds_and_cs(2)
+#H = get_1body_ham(h, cds, cs)
+#H_N = get_N_subspace(H, Nop, N=2)
+#N_N = get_N_subspace(Nop, Nop, N=2)
+#print np.around(H_N,1)
+#print np.around(N_N,1)
+# three basis functions
+h = -2.0 * (np.eye(3,3,1) + np.eye(3,3,-1))
+cds, cs, Nop = get_cds_and_cs(3)
+print cds[0][0].shape
+H = get_1body_ham(h, cds, cs)
+H_N = get_N_subspace(H, Nop, N=3)
+N_N = get_N_subspace(Nop, Nop, N=3)
+print np.around(H_N,1)
+print np.around(N_N,1)
+

Added: trunk/kgf/exact_diagonalization/exactd3.py
===================================================================
--- trunk/kgf/exact_diagonalization/exactd3.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/exactd3.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,340 @@
+import numpy as np
+from kgf.Tools import dots
+from numpy import ascontiguousarray as asc
+from numpy import linalg as la
+#import pylab as pl
+
+# one single particle state
+class ExactDiag:
+    def __init__(self, nstates=1):
+        self.P1 = np.identity(4,int)
+        self.S1z = np.zeros((4,4),int)
+        self.c1da = np.zeros((4,4),int)
+        self.c1db = np.zeros((4,4),int)
+        self.c1da[1,0] = 1
+        self.c1da[3,2] = 1
+        self.c1db[2,0] = 1
+        self.c1db[3,1] = -1
+        self.c1a = asc(self.c1da.T)
+        self.c1b = asc(self.c1db.T)
+        self.P1[1,1] = -1
+        self.P1[2,2] = -1
+        self.S1z[1,1] = 1
+        self.S1z[2,2] = -1
+        self.N1 = np.dot(self.c1da, self.c1a) + np.dot(self.c1db, self.c1b)
+        self.I1 = np.identity(4, int)
+
+        self.N = self.N1.copy()
+        self.Sz = self.S1z.copy()
+        self.P = self.P1.copy()
+        self.cd_si = [[self.c1da.copy()], [self.c1db.copy()]]
+        self.c_si = [[self.c1a.copy], [self.c1b.copy()]]
+        #self.cda = [self.c1da.copy()] # c^dagger_i\uparrow
+        #self.cdb = [self.c1db.copy()] # c^dagger_i\downarrow
+        #self.ca = [self.c1a.copy()]
+        #self.cb = [self.c1b.copy()]
+        self.n = 1
+        self.onebody_Nwvi = {}
+
+        for i in range(nstates-1):
+            self.add_new_state()
+
+    def set_onebody_ham(self, h):
+        self.h = self.get_onebody_ham(h)
+
+    def diagonalize_onebody_ham(self, Npart):
+        indices = self.get_Npart_subspace(Npart)
+        h = self.h.take(indices, axis=0).take(indices, axis=1)
+        w, v = la.eig(h)
+        sortlist = w.argsort()
+        w = w[sortlist]
+        v = v.take(sortlist, axis=1)
+        self.onebody_Nwvi[Npart] = {'v':v, 'w':w, 'indices':indices}
+   
+    def get_overlap(self,N1, N2, M1, M2):
+        Nwvi = self.onebody_Nwvi
+        if N1 not in Nwvi.keys():
+            self.diagonalize_onebody_ham(N1)
+        if N2 not in Nwvi.keys():
+            self.diagonalize_onebody_ham(N2)
+
+        indices1 = Nwvi[N1]['indices'] 
+        indices2 = Nwvi[N2]['indices']
+
+    def get_homo_energy(self, Nneutral):
+        Nwvi = self.onebody_Nwvi 
+        E = Nwvi[Nneutral]['w'][0]
+        Em1 = Nwvi[Nneutral-1]['w'][0]
+        return E - Em1
+
+    def get_lumo_energy(self, Nneutral):
+        Nwvi = self.onebody_Nwvi 
+        E = Nwvi[Nneutral]['w'][0]
+        Ep1 = Nwvi[Nneutral+1]['w'][0]
+        return Ep1 - E
+
+    def get_homo_lumo_gap(self, Nneutral):
+        el = self.get_lumo_energy(Nneutral)
+        eh = self.get_homo_energy(Nneutral)
+        return el - eh
+
+    def get_number_of_states(self):
+        return self.n
+
+    def get_size(self):
+        return self.cd_si[0][0].shape[0]
+
+    def add_new_state(self):
+        n = self.get_number_of_states()
+        for i in range(n):
+            for s in range(2):
+                self.cd_si[s][i] = np.kron(self.cd_si[s][i], self.I1)
+                self.c_si[s][i] = asc(self.cd_si[s][i].T)
+
+        cdanew = np.kron(self.P, self.c1da)
+        cdbnew = np.kron(self.P, self.c1db)
+        canew = asc(cdanew.T)
+        cbnew = asc(cdbnew.T)
+        
+        self.cd_si[0].append(cdanew)
+        self.cd_si[1].append(cdbnew)
+        self.c_si[0].append(canew)
+        self.c_si[1].append(cbnew)
+
+        self.P = np.kron(self.P, self.P1)
+        self.N = (np.kron(self.N, self.I1) + 
+                  np.kron(np.identity(len(self.N),int), self.N1))
+        
+        self.Sz = (np.kron(self.Sz, self.I1) + 
+                   np.kron(np.identity(len(self.Sz),int), self.S1z))
+        self.n += 1
+
+    def get_N(self, debug=False):
+        if not debug:
+            return self.N
+        else:
+            n = self.get_number_of_states()
+            size = self.get_size()
+            N = np.zeros((size,size), int)
+            for i in range(n):
+                N += np.dot(self.cdas[i], self.cas[i])    
+                N += np.dot(self.cdbs[i], self.cbs[i])
+
+            return N
+
+    def get_Npart_subspace(self, N):
+        return [i for i, n in enumerate(self.N.diagonal()) if n==N]
+
+    def get_Sz_subspace(self, Sz):
+        sz = int(2 * Sz)
+        return [i for i, s in enumerate(self.Sz.diagonal()) if s==sz]
+
+    def get_Npart_Sz_subspace(self, N, Sz):
+        sz = int(2 * Sz)
+        indices = self.get_Npart_subspace(N)
+        Sz_sub = self.Sz.diagonal()[indices]
+        return [i for i, sz_sub in zip(indices, Sz_sub) if sz_sub==sz]
+
+    def get_onebody_ham(self, h, N=None, Sz=None):
+        n = h.shape[0]
+        assert n==self.get_number_of_states()
+        H = self.zeros(h.dtype)
+        cd_si, c_si = self.cd_si, self.c_si
+        for i in range(n):
+            for j in range(n):
+                for s in range(2):
+                    H += h[i,j] * np.dot(cd_si[s][i], c_si[s][j])
+       
+        if N is not None:
+            indices = self.get_Npart_subspace(N)
+            H = asc(H.take(indices, axis=0).take(indices,axis=1))
+        return H
+
+    def get_ni(self, i, spin=None):
+        """
+            out: n_i = c^\dagger_{i,spin} c^_{i,spin}
+            spin : int {0, 1}, None
+        """
+        n_i = self.zeros(int)
+        if spin is None:
+            spin = [0, 1]
+        else:
+            spin = [spin]
+       
+        for s in spin:
+            n_i += dots(self.cd_si[s][i], self.c_si[s][i]) 
+    
+        return n_i
+
+    def get_N2(self):
+        N = self.zeros(int)
+        for i in range(self.get_number_of_states()):
+            N += self.get_ni(i)
+        
+        return N
+           
+    def zeros(self, dtype=float):
+        n = self.get_number_of_states()
+        return np.zeros((4**n, 4**n), dtype)
+    
+    def get_Sz_S2_S_values(self, Npart):
+        data = self.onebody_Nwvi[Npart]
+        indices, v = data['indices'], data['v']
+        Sz = self.get_Sz()
+        Sz = Sz.take(indices, axis=0).take(indices, axis=1)
+        Sz = dots(v.T, Sz, v).diagonal()
+        S2 = self.get_S2()
+        S2 = S2.take(indices, axis=0).take(indices, axis=1)
+        S2 = dots(v.T, S2, v).diagonal()
+        S = -0.5 + np.sqrt(0.25+S2)
+        return Sz, S2, S 
+
+    def get_Sz(self, power=None):
+        n = self.get_number_of_states()
+        Sz = self.zeros()
+        for i in range(n):
+            Sz += self.get_ni(i, 0) - self.get_ni(i, 1)        
+        
+        Sz *= 0.5
+        if power==2:
+            return dots(Sz, Sz)
+        else:
+            return Sz
+ 
+    def get_Sx(self, power=None):
+        n = self.get_number_of_states()
+        Sx = self.zeros()
+        cd, c = self.cd_si, self.c_si
+        for i in range(n):
+            Sx += dots(cd[0][i], c[1][i]) + dots(cd[1][i], c[0][i])
+            #Sx += dots(self.cda[i], self.cb[i]) + dots(self.cdb[i], self.ca[i])
+
+        Sx *= 0.5
+        if power==2:
+            return dots(Sx,Sx)
+        else:    
+            return Sx
+
+    def get_Sy(self, power=None):
+        n = self.get_number_of_states()
+        Sy = self.zeros()
+        cd, c = self.cd_si, self.c_si
+        for i in range(n):
+            Sy += dots(cd[1][i], c[0][i]) - dots(cd[0][i], c[1][i])
+            #Sy += dots(self.cdb[i], self.ca[i]) - dots(self.cda[i], self.cb[i])
+      
+        if power==2:
+            return -0.25 * dots(Sy,Sy) 
+        else:
+            return Sy * 0.5j
+  
+    def get_S2(self):
+        S2 = self.get_Sx(2) + self.get_Sy(2) + self.get_Sz(2)
+        return S2
+
+    def print_occupation_numbers(self, state):
+        """
+            print occupations of single particle states 
+            prints: i,  spin, <state|n_{i\sigma}|state> 
+        """
+        if type(state)==int:
+            state = np.eye(1, self.get_size(), state).T
+
+        n = self.get_number_of_states()
+        cd, c = self.cd_si, self.c_si
+        print '-----------'
+        print 'i  s <n_is>'
+        print '-----------'
+        for i in range(n):
+            for s1 in range(2):
+                nop = dots(cd[s1][i], c[s1][i])
+                nval = dots(state.T, nop, state)
+                if abs(nval)>1.0e-15:
+                    print '%.2i %i %.3f' % (i, s1, nval)
+   
+    def get_hole_overlap(self, N1, N2, M1, M2, i, spin):
+        """ 
+             h_M1(i,spin) = <psi^{N1}_M1 | c_{i,spin} | psi^{N2}_M2 >
+             usually N2 = N_neutral
+                     M2 = 0 
+                     N1 = N2-1
+        """          
+        assert N1 and N2 in self.onebody_Nwvi.keys()
+        data1 = self.onebody_Nwvi[N1]
+        data2 = self.onebody_Nwvi[N2]
+        c = self.c_si[spin][i]
+        c = c.take(data1['indices'], axis=0).take(data2['indices'], axis=1)
+        psiN1_M1 = data1['v'][:,M1]
+        psiN2_M2 = data2['v'][:,M2]
+        h_M1 = dots(psiN1_M1.T, c, psiN2_M2)
+        return h_M1
+
+    def get_hubbard_ham(self, U):
+        n = self.get_number_of_states()
+        assert n==U.shape[0]
+        V = self.zeros()
+        for i in range(n):
+            V += U[i,i] * np.dot(self.get_ni(i,0), self.get_ni(i,1))
+        
+        return V
+            
+    def get_2body_ham(self, v):
+            """       ___
+                    1 \         +   +
+                V = -  \    V  c   c   c    c
+                    2 /      ij is  js' js'  is
+                     /___     
+                        ij,s,s'
+            """
+            n = h.shape[0]
+            V = np.zeros((4**n, 4**n), v.dtype)
+            for i in range(n):
+                for j in range(n):
+                    # ia ja
+                    cdia, cia = self.cda[i], self.ca[i] # c^+_{i+}, c_{i+}
+                    cdja, cja = self.cda[j], self.ca[j] # c^_{j+}, c_{j+}
+                    operator = np.dot(np.dot(np.dot(cdia, cdja), cja), cia)
+                    V += (0.5 * v[i,j]) * operator
+                    # ia jb
+                    cdia, cia = self.cda[i], self.ca[i] # c^+_{i+}, c_{i+}
+                    cdja, cja = self.cda[j], self.ca[j] # c^_{j+}, c_{j+}
+                    operator = np.dot(np.dot(np.dot(cdia, cdja), cja), cia)
+                    V += (0.5 * v[i,j]) * operator
+                    #b a
+               
+
+#h = -2.0*(np.eye(2,2,1) + np.eye(2,2,-1))
+
+if __name__ == '__main__':
+    calc = ExactDiag(2)
+    h = np.zeros((2, 2))
+    U = np.identity(2) 
+    t = -1.0
+    u = 2 * abs(t)
+    h[0, 1] = t
+    h[1, 0] = t
+    U *= u
+    H = calc.get_onebody_ham(h)
+    V = calc.get_hubbard_ham(U)
+    H2 = H + V
+    indices = calc.get_Npart_subspace(2)
+    H2sub = asc(H2.take(indices,axis=0).take(indices,axis=1))
+    Nsub = asc(calc.N.take(indices,axis=0).take(indices,axis=1))
+    S2sub = asc(calc.get_S2().take(indices,axis=0).take(indices,axis=1))
+    Szsub = asc(calc.get_Sz().take(indices,axis=0).take(indices,axis=1))
+
+    eigs,vecs = la.eigh(H2sub)
+    print np.around(eigs,3)
+    print np.around(dots(vecs.T,Nsub,vecs),3).diagonal()
+    print np.around(dots(vecs.T,Szsub,vecs),3).diagonal()
+    print np.around(dots(vecs.T,S2sub,vecs),3).diagonal()
+    #print H
+#    eigs, vecs = np.linalg.eigh(H)
+#    print np.around(eigs,1)
+#    print np.around(vecs[:,0],2)
+    # configurations
+#    conf = np.zeros(16)
+#    conf[0] = 1.0
+    # |b+,b-> 
+    # |b,s> = 1 / sqrt(2) (|1,s> + |2,s>)
+    # 

Added: trunk/kgf/exact_diagonalization/exactdsparse.py
===================================================================
--- trunk/kgf/exact_diagonalization/exactdsparse.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/exactdsparse.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,239 @@
+from sparsetools import get_subspace
+from scipy import sparse as sp
+from scipy.sparse import linalg as spla
+import numpy as np
+from numpy import linalg as npla
+from numpy import int8
+
+class ExactDiag:
+    def __init__(self, nsites):
+        self.P1 = sp.csr_matrix([[1, 0, 0, 0],
+                                 [0,-1, 0, 0],
+                                 [0, 0,-1, 0],
+                                 [0, 0, 0, 1]], dtype=int8)
+
+        self.S1z = sp.csr_matrix([[0, 0, 0, 0],
+                                  [0, 1, 0, 0],
+                                  [0, 0,-1, 0],
+                                  [0, 0, 0, 0]], dtype=int8)
+
+        self.c1da = sp.csr_matrix([[0, 0, 0, 0],
+                                   [1, 0, 0, 0],
+                                   [0, 0, 0, 0],
+                                   [0, 0, 1, 0]], dtype=int8)
+
+        self.c1db = sp.csr_matrix([[0, 0, 0, 0],
+                                   [0, 0, 0, 0],
+                                   [1, 0, 0, 0],
+                                   [0,-1, 0, 0]], dtype=int8)
+
+        self.c1a = self.c1da.T
+        self.c1b = self.c1db.T
+        self.I1 = sp.identity(4, dtype=int8, format='csr')
+
+        self.P = self.P1.copy()
+        self.cd_si = [[self.c1da.copy()], [self.c1db.copy()]]
+        self.c_si = [[self.c1a.copy()], [self.c1b.copy()]]
+        self.data_Nwvi = {}
+
+        self.nsites = 1
+        self.shape = (4**self.nsites, 4**self.nsites)
+        for i in range(nsites-1):
+            self.add_new_state()
+          
+    def add_new_state(self):
+        for i in range(self.nsites):
+            for s in range(2):
+                self.cd_si[s][i] = sp.kron(self.cd_si[s][i], self.I1, 'csr')
+                self.c_si[s][i] = self.cd_si[s][i].T
+
+        cdanew = sp.kron(self.P, self.c1da)#, 'csr')
+        cdbnew = sp.kron(self.P, self.c1db)#, 'csr')
+        
+        self.cd_si[0].append(cdanew)
+        self.cd_si[1].append(cdbnew)
+        self.c_si[0].append(cdanew.T)
+        self.c_si[1].append(cdbnew.T)
+
+        self.P = sp.kron(self.P, self.P1, 'csr')
+        self.nsites += 1
+        self.shape = (4**self.nsites, 4**self.nsites)
+
+    def get_N(self):
+        N = sp.csr_matrix(self.shape, dtype=int8)
+        cd, c = self.cd_si, self.c_si
+        for i in range(self.nsites):
+            for s in range(2):
+                N = N + cd[s][i].dot(c[s][i])
+        return N
+   
+    def get_Sx(self):
+        Sx = sp.csr_matrix(self.shape, dtype=int8)
+        cd, c = self.cd_si, self.c_si
+        for i in range(self.nsites):
+            Sx = Sx + cd[0][i].dot(c[1][i]) + cd[1][i].dot(c[0][i])
+        return Sx, 0.5
+
+    def get_Sy(self):
+        Sy = sp.csr_matrix(self.shape, dtype=int8)
+        cd, c = self.cd_si, self.c_si
+        for i in range(self.nsites):
+            Sy = Sy + cd[0][i].dot(c[1][i]) - cd[1][i].dot(c[0][i])
+        return Sy, -0.5j
+
+    def get_Sz(self):
+        Sz = sp.csr_matrix(self.shape, dtype=int8)
+        cd, c = self.cd_si, self.c_si
+        for i in range(self.nsites):
+            Sz = Sz + cd[0][i].dot(c[0][i]) - cd[1][i].dot(c[1][i])
+        return Sz, 0.5
+      
+    def get_S2(self):
+        Sx, px = self.get_Sx()
+        Sy, py = self.get_Sy()
+        Sz, pz = self.get_Sz()
+        S2 = Sx.dot(Sx) - Sy.dot(Sy) + Sz.dot(Sz)
+        return S2, 0.25
+    
+    def get_subspace_indices(self, npart=None, sz=None, s=None):
+        """
+            indices for subspace  
+        """
+        if npart is not None:
+            ndia = self.get_N().diagonal()
+            indices = [i for i, n in enumerate(ndia) if n==npart]
+        else:
+            indices = range(self.shape[0])
+        
+        if sz is not None:
+            szi = int(2 * sz)
+            sz2, pref = self.get_Sz()
+            indices = [i for i in indices if sz2[i,i]==szi]
+       
+        if s is not None:
+            assert 0, 'not implemented'
+            s2i = int(4*s*(s+1))
+            #s22, pref = self.get_S2()
+            #indices = [i for i in indices if s22[i,i]==s2i]
+
+        return indices
+    
+    def get_subspace(self, A, indices):
+        return get_subspace(A, indices, 'csr')
+
+    def get_onebody_ham(self, h):
+        cd, c = self.cd_si, self.c_si
+        h2 = sp.lil_matrix(self.shape, dtype=h.dtype)
+        for i in range(self.nsites):
+            for j in range(self.nsites):
+                for s in range(2):
+                    h2 = h2 + h[i, j] * cd[s][i].dot(c[s][j])
+        return h2
+
+    def get_twobody_hubbard_ham(self, v):
+        cd, c = self.cd_si, self.c_si
+        v2 = sp.lil_matrix(self.shape, dtype=v.dtype)
+        for i in range(self.nsites):
+            nia = cd[0][i].dot(c[0][i])
+            nib = cd[1][i].dot(c[1][i])
+            v2 = v2 + v[i, i] * nia.dot(nib)
+        return v2
+
+    def get_twobody_PPP_ham(self, v):
+        cd, c = self.cd_si, self.c_si
+        v2 = sp.lil_matrix(self.shape, dtype=v.dtype)
+        for i in range(self.nsites):
+            for j in range(self.nsites):
+                if i==j: 
+                    nia = cd[0][i].dot(c[0][i])
+                    nib = cd[1][i].dot(c[1][i])
+                    v2 = v2 + v[i, i] * nia.dot(nib)
+                else:
+                    ni = cd[0][i].dot(c[0][i]) + cd[1][i].dot(c[1][i])
+                    nj = cd[0][j].dot(c[0][j]) + cd[1][j].dot(c[1][j])
+                    v2 = v2 + (0.5 * v[i, j]) * ni.dot(nj)
+        return v2
+
+    def get_expectation_value(self, R, vec, diag_only=False):
+        if type(vec)==np.ndarray:
+            vec = sp.csr_matrix(vec)
+        out = vec.T.dot(R).dot(vec)
+        if diag_only:
+            out = out.diagonal()
+        
+        return out
+
+    def get_lowest_energy(self, H):
+        val, vec = spla.eigsh(H)
+        return val[0]
+
+    def get_expectation_values(self, H, indices):
+        h_nn = self.get_subspace(H, indices)
+        if 0:
+            val, vec = spla.eigsh(h_nn, k=h_nn.shape[0]-1)
+        else:
+            val, vec = npla.eig(h_nn.toarray())
+            perm = val.argsort()
+            val = val[perm]
+            vec = vec.take(perm, axis=1)
+        Sz_nn = self.get_subspace(self.get_Sz()[0], indices)
+        S2_nn = self.get_subspace(self.get_S2()[0], indices)
+        n_nn = self.get_subspace(self.get_N(), indices)
+        n_n = self.get_expectation_value(n_nn, vec, True)
+        e_n = self.get_expectation_value(h_nn, vec, True)
+        sz_n = self.get_expectation_value(Sz_nn, vec, True) * 0.5
+        s2_n = self.get_expectation_value(S2_nn, vec, True) * 0.25
+        s_n = -0.5 + np.sqrt(s2_n + 0.25)
+        return e_n, sz_n, s_n, s2_n, n_n
+
+    def print_expectation_values(self, H, indices, nmax=5):
+        e_n, sz_n, s_n, s2_n, n_n = self.get_expectation_values(H, indices)
+        print '--------------------------------'
+        print '  E     |  Sz   |  S    |   N'
+        print '--------------------------------'
+        for e, sz, s2, n in zip(e_n[:nmax], 
+                                sz_n[:nmax], s_n[:nmax], n_n[:nmax]):
+            es = '%.3f |' % e
+            szs = '%.2f |' % sz
+            s2s = '%.2f |' % s2
+            nvs = '%.2f' % n
+            print es.rjust(9), szs.rjust(7), s2s.rjust(7), nvs.rjust(5)
+
+    def print_occupation_numbers(self, state):
+        """
+            print occupations of single particle states 
+            prints: i,  spin, <state|n_{i\sigma}|state> 
+        """
+        if type(state)==int:
+            state = sp.eye(1, self.shape[0], state, dtype=int8).T
+
+        cd, c = self.cd_si, self.c_si
+        print '-----------'
+        print 'i  s <n_is>'
+        print '-----------'
+        for i in range(self.nsites):
+            for s in range(2):
+                nop = cd[s][i].dot(c[s][i])
+                n = state.T.dot(nop).dot(state)[0,0]
+                if abs(n)>1.0e-15:
+                    print '%.2i %i %.3f' % (i, s, n)
+
+if __name__=='__main__':
+    calc = ExactDiag(2)
+    Sz, pz = calc.get_Sz()
+    S2, p2 = calc.get_S2()
+
+    h1 = np.array([[ 0.0,-1.0],
+                   [-1.0, 0.0]])
+    V1 = np.identity(2) * 2.0
+    V2 = calc.get_twobody_hubbard_ham(V1)
+    h2 = calc.get_onebody_ham(h1)
+    indices = calc.get_subspace_indices(npart=2)# sz=0.0)
+    V2t = calc.get_subspace(V2, indices)
+    h2t = calc.get_subspace(h2, indices)
+    Szt = calc.get_subspace(Sz, indices)
+    S2t = calc.get_subspace(S2, indices)
+    H = h2t + V2t
+    val,vec = npla.eigh(H.toarray())
+    Szb = calc.get_expectation_value(Szt, vec)
+    S2b = calc.get_expectation_value(S2t, vec)

Added: trunk/kgf/exact_diagonalization/sparsetools.py
===================================================================
--- trunk/kgf/exact_diagonalization/sparsetools.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/sparsetools.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,25 @@
+from scipy import sparse as sp
+
+def get_subspace(A, indices, form='csr'):
+    n = len(indices)
+    At = A.tolil()
+    B = sp.lil_matrix((n,n), dtype=A.dtype)
+    for i in range(n):
+        for j in range(i, n):
+            B[i, j] = At[indices[i], indices[j]]
+
+    diag = B.diagonal()
+    B = (B + B.T).tolil()
+    B.setdiag(diag)
+    return B.asformat(form)
+
+def get_subspace2(A, indices, form='csr'):
+    n = len(indices)
+    At = A.tolil()
+    B = sp.lil_matrix((n,n), dtype=A.dtype)
+    for i, i1 in enumerate(indices):
+        for j, j1 in enumerate(indices):
+            B[i, j] = At[i1, j1] 
+
+    return B.asformat(form)
+

Added: trunk/kgf/exact_diagonalization/test.py
===================================================================
--- trunk/kgf/exact_diagonalization/test.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/test.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,14 @@
+from kgf.Tools import dots
+import numpy as np
+from numpy import linalg as la
+from exactd3 import ExactDiag
+
+calc = ExactDiag(2)
+indices = calc.get_Npart_subspace(2)
+Sz = calc.get_Sz().take(indices,axis=0).take(indices, axis=1)
+S2 = calc.get_S2().take(indices,axis=0).take(indices, axis=1)
+
+w2, v2 = la.eigh(S2)
+x = dots(v2.conj().T, Sz, v2)
+print 'Sz:', x.diagonal()
+print 'S2:', w2

Added: trunk/kgf/exact_diagonalization/test_hubbard.py
===================================================================
--- trunk/kgf/exact_diagonalization/test_hubbard.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/test_hubbard.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,69 @@
+from kgf.Tools import dots
+import numpy as np
+from numpy import around, ascontiguousarray as asc
+from numpy import linalg as la
+from exactd3 import ExactDiag
+#import pylab as pl
+
+t = -1.0
+u = 2 * abs(t)
+h1b = np.zeros((2, 2))
+#h1b[0, 1] = t
+#h1b[1, 0] = t
+h1b[0,0] = -1
+h1b[1,1] = 1
+U = np.identity(2) * u
+
+calc = ExactDiag(2)
+h = calc.get_onebody_ham(h1b)
+V = calc.get_hubbard_ham(U)
+Ssqr = calc.get_S2()
+Sz = calc.get_Sz()
+N = calc.get_N()
+Htot = h + V
+indices = calc.get_Npart_subspace(2)
+h2 = h.take(indices,axis=0).take(indices,axis=1)
+V2 = V.take(indices,axis=0).take(indices,axis=1)
+Htot2 = h2 + V2
+N2 = N.take(indices,axis=0).take(indices,axis=1)
+Ssqr2 = Ssqr.take(indices,axis=0).take(indices,axis=1)
+Sz2 = Sz.take(indices,axis=0).take(indices,axis=1)
+
+w, v = la.eig(h)
+Szval = dots(v.T.conj(), Sz, v)
+Ssqrval = dots(v.T.conj(), Ssqr, v)
+
+print 'En:\n', np.around(w,2)
+print 'Sz:\n', np.around(Szval.diagonal(),2)
+print 'S2:\n', np.around(Ssqrval.diagonal(),2)
+
+calc.set_onebody_ham(h1b)
+calc.diagonalize_onebody_ham(1)
+calc.diagonalize_onebody_ham(2)
+calc.diagonalize_onebody_ham(3)
+
+calc.print_occupation_numbers(2)
+
+if 0:
+    gs = v[:, w.argmin()]
+    cb = (1.0 / np.sqrt(2)) * (calc.ca[0] + calc.ca[1])
+    cgs = dots(cb,gs)
+    dots(cgs.T,Sz,cgs)
+
+    w2, v2 = la.eig(h2)
+    Szval2 = dots(v2.T.conj(), Sz2, v2)
+    Ssqrval2 = dots(v2.T.conj(), Ssqr2, v2)
+
+    print 'En2:\n', np.around(w2,2)
+    print 'Sz2:\n', np.around(Szval2.diagonal(),2)
+    print 'S22:\n', np.around(Ssqrval2.diagonal(),2)
+
+
+
+
+if 0:
+    pl.imshow(Szval, interpolation='nearest')
+    pl.show()
+    pl.imshow(Ssqrval, interpolation='nearest')
+    pl.show()
+

Added: trunk/kgf/exact_diagonalization/test_hubbard_sparse.py
===================================================================
--- trunk/kgf/exact_diagonalization/test_hubbard_sparse.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/test_hubbard_sparse.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,30 @@
+import numpy as np
+from numpy import linalg as npla
+from scipy import sparse as sp
+from scipy.sparse import linalg as spla
+from exactdsparse import ExactDiag
+
+calc = ExactDiag(2)
+Sz, pz = calc.get_Sz()
+S2, p2 = calc.get_S2()
+N = calc.get_N()
+t = -5.0
+u = 10.0
+h1 = np.array([[ 0.0+1.0e-10, t],
+               [   t, 0.0]])               
+V1 = np.identity(2) * u
+h_ion = np.zeros((2,2))
+for i in range(2):
+    h_ion[i, i] += -0.5 * V1[i,i]
+    for j in range(2):
+        if j!=i:
+            h_ion[i,j] += 1.0 * V1[i,j]
+
+
+V2 = calc.get_twobody_hubbard_ham(V1)
+h2 = calc.get_onebody_ham(h1)
+H = h2 + V2
+for npart in [1,2,3]:
+    indices = calc.get_subspace_indices(npart=npart)
+    calc.print_expectation_values(H, indices)
+

Added: trunk/kgf/exact_diagonalization/test_ppp_sparse.py
===================================================================
--- trunk/kgf/exact_diagonalization/test_ppp_sparse.py	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/test_ppp_sparse.py	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,32 @@
+import numpy as np
+from numpy import linalg as npla
+from scipy import sparse as sp
+from scipy.sparse import linalg as spla
+from exactdsparse import ExactDiag
+
+calc = ExactDiag(2)
+Sz, pz = calc.get_Sz()
+S2, p2 = calc.get_S2()
+N = calc.get_N()
+t = -8.0
+u = 10.0
+e0 = 0.0
+h1 = np.array([[ e0+1.0e-5, t],
+               [   t, e0]])               
+
+V1 = np.identity(2) * u
+V1[0,1] = V1[1,0] = 2*u/3.
+h_ion = np.zeros((2,2))
+for i in range(2):
+    h_ion[i, i] += -0.5 * V1[i,i]
+    for j in range(2):
+        if j!=i:
+            h_ion[i,i] += -1.0 * V1[i,j]
+
+V2 = calc.get_twobody_hubbard_ham(V1)
+h2 = calc.get_onebody_ham(h1 + h_ion)
+H = h2 + V2
+for npart in [1,2,3]:
+    indices = calc.get_subspace_indices(npart=npart)
+    calc.print_expectation_values(H, indices)
+

Added: trunk/kgf/exact_diagonalization/tmp.dat
===================================================================
--- trunk/kgf/exact_diagonalization/tmp.dat	                        (rev 0)
+++ trunk/kgf/exact_diagonalization/tmp.dat	2012-11-16 14:10:50 UTC (rev 268)
@@ -0,0 +1,4 @@
+200 0.7
+400 2.2
+800 7.9
+1600 30.9



More information about the KeldyshGF-svncheckins mailing list