[KeldyshGF-svncheckins] [svn commit /var/www/svn/KeldyshGF] r270 - trunk/kgf/exact_diagonalization

svn at fysik.dtu.dk svn at fysik.dtu.dk
Fri Nov 16 15:13:43 CET 2012


Author: strange
Date: 2012-11-16 15:13:43 +0100 (Fri, 16 Nov 2012)
New Revision: 270

Removed:
   trunk/kgf/exact_diagonalization/exactd2.py
   trunk/kgf/exact_diagonalization/exactd3.py
   trunk/kgf/exact_diagonalization/test_hubbard.py
   trunk/kgf/exact_diagonalization/test_hubbard_sparse.py
   trunk/kgf/exact_diagonalization/tmp.dat
Log:
arghh...

Deleted: trunk/kgf/exact_diagonalization/exactd2.py
===================================================================
--- trunk/kgf/exact_diagonalization/exactd2.py	2012-11-16 14:13:06 UTC (rev 269)
+++ trunk/kgf/exact_diagonalization/exactd2.py	2012-11-16 14:13:43 UTC (rev 270)
@@ -1,81 +0,0 @@
-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)
-

Deleted: trunk/kgf/exact_diagonalization/exactd3.py
===================================================================
--- trunk/kgf/exact_diagonalization/exactd3.py	2012-11-16 14:13:06 UTC (rev 269)
+++ trunk/kgf/exact_diagonalization/exactd3.py	2012-11-16 14:13:43 UTC (rev 270)
@@ -1,340 +0,0 @@
-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>)
-    # 

Deleted: trunk/kgf/exact_diagonalization/test_hubbard.py
===================================================================
--- trunk/kgf/exact_diagonalization/test_hubbard.py	2012-11-16 14:13:06 UTC (rev 269)
+++ trunk/kgf/exact_diagonalization/test_hubbard.py	2012-11-16 14:13:43 UTC (rev 270)
@@ -1,69 +0,0 @@
-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()
-

Deleted: trunk/kgf/exact_diagonalization/test_hubbard_sparse.py
===================================================================
--- trunk/kgf/exact_diagonalization/test_hubbard_sparse.py	2012-11-16 14:13:06 UTC (rev 269)
+++ trunk/kgf/exact_diagonalization/test_hubbard_sparse.py	2012-11-16 14:13:43 UTC (rev 270)
@@ -1,30 +0,0 @@
-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)
-

Deleted: trunk/kgf/exact_diagonalization/tmp.dat
===================================================================
--- trunk/kgf/exact_diagonalization/tmp.dat	2012-11-16 14:13:06 UTC (rev 269)
+++ trunk/kgf/exact_diagonalization/tmp.dat	2012-11-16 14:13:43 UTC (rev 270)
@@ -1,4 +0,0 @@
-200 0.7
-400 2.2
-800 7.9
-1600 30.9



More information about the KeldyshGF-svncheckins mailing list