Merge pull request #92 from freude/self_energies

Self energies

examples/surf_greens_function.py

@@ -1,7 +1,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import nanonet.tb as tb
-from nanonet.negf.greens_functions import simple_iterative_greens_function, surface_greens_function
+from nanonet.negf.greens_functions import simple_iterative_greens_function, sancho_rubio_iterative_greens_function, surface_greens_function 
 
 
 def main(surf_greens_fun):
@@ -30,13 +30,13 @@ def main(surf_greens_fun):
     sgf_r = []
 
     for E in energy:
-        sf = surf_greens_fun(E, h_l, h_0, h_r)
+        sf = surf_greens_fun(E, h_l, h_0, h_r , damp=0.001j)
         if isinstance(sf, tuple):
             L = sf[0]
             R = sf[1]
         else:
             L = sf
-            R = surf_greens_fun(E, h_r, h_0, h_l)
+            R = surf_greens_fun(E, h_r, h_0, h_l , damp=0.001j)
 
         sgf_l.append(L)
         sgf_r.append(R)
@@ -108,13 +108,13 @@ def main1(surf_greens_fun):
     sgf_r = []
 
     for E in energy:
-        sf = surf_greens_fun(E, h_l, h_0, h_r)
+        sf = surf_greens_fun(E, h_l, h_0, h_r, damp=0.001j)
         if isinstance(sf, tuple):
             L = sf[0]
             R = sf[1]
         else:
             L = sf
-            R = surf_greens_fun(E, h_r, h_0, h_l)
+            R = surf_greens_fun(E, h_r, h_0, h_l, damp=0.001j)
 
         sgf_l.append(L)
         sgf_r.append(R)
@@ -153,6 +153,8 @@ def main1(surf_greens_fun):
 
     main(surf_greens_fun=surface_greens_function)
     main(surf_greens_fun=simple_iterative_greens_function)
+    main(surf_greens_fun=sancho_rubio_iterative_greens_function)
     main1(surf_greens_fun=surface_greens_function)
     main1(surf_greens_fun=simple_iterative_greens_function)
+    main1(surf_greens_fun=sancho_rubio_iterative_greens_function)
 

nanonet/negf/greens_functions.py

@@ -286,21 +286,21 @@ def simple_iterative_greens_function(E, h_l, h_0, h_r, **kwargs):
             s_r : numpy array (Default value = zero matrix)
                   right block of three-block-diagonal overlap matrix
 
-            alpha : float (Default value = 0.5)
-                    mixing parameter (0 < alpha < 1) for green's function
+          alpha : float (Default value = 0.5)
+                  mixing parameter (0 < alpha < 1) for green's function
 
-            posinf : float (Default value = 1e-6)
-                     positive infinitesimal a.k.a broadening parameter
-                     that shifts solution off real axis
+           damp : float (Default value = 1e-6)
+                  positive infinitesimal a.k.a broadening parameter
+                  that shifts solution off real axis
 
-            initialguess : numpy array (Default value = zero matrix)
-                     initial guess of the self-energy, allows restarting of algorithm
-                     to gain increased precision, or if the user has a good idea of
-                     what the self-energy should be, can boost convergence
+   initialguess : numpy array (Default value = zero matrix)
+                  initial guess of the self-energy, allows restarting of algorithm
+                  to gain increased precision, or if the user has a good idea of
+                  what the self-energy should be, can boost convergence
 
-            nconv : float (Default value = 1e-10)
-                    convergence requirement for iteration using the 2-norm,
-                    also known as the largest singular value.
+          nconv : float (Default value = 1e-10)
+                  convergence requirement for iteration using the 2-norm,
+                  also known as the largest singular value.
 
     Returns
     -------
@@ -313,11 +313,11 @@ def simple_iterative_greens_function(E, h_l, h_0, h_r, **kwargs):
         and convergence criterion in order to compute the retarded self-energy.
 
         Note that due to its singular nature, computing the self-energy of a
-        matrix with posinf=0 is numerically undefined when starting from scratch,
+        matrix with damp=0 is numerically undefined when starting from scratch,
         however, if the user desires, they may compute the solution at some small
-        posinf, say 1e-8, then use that as an input guess for a posinf=0 input,
+        damp, say 1e-8, then use that as an input guess for a damp=0 input,
         the imaginary part of the initial guess shifts the solution off the real
-        line and allows for convergence. We do not suggest starting a posinf=0
+        line and allows for convergence. We do not suggest starting a damp=0
         iteration from an arbitrary guess due to instability in convergence.
 
     """
@@ -326,26 +326,24 @@ def simple_iterative_greens_function(E, h_l, h_0, h_r, **kwargs):
     s_0 = kwargs.get('s_0', np.identity(h_0.shape[0]))
     s_r = kwargs.get('s_r', np.zeros(h_0.shape[0]))
     alpha = kwargs.get('alpha', 0.5)
-    posinf = kwargs.get('posinf', 1e-6)
+    damp = kwargs.get('damp', 1e-6)
     initialguess = kwargs.get('initialguess', np.zeros(h_0.shape[0]))
     nconv = kwargs.get('nconv', 1e-10)
 
-    #####This is where the function would start if working properly:
-
     # Check inputs:
     # If alpha is in an unphysical range, fix it to 0.5:
     if (alpha < 0) or (alpha > 1):
         alpha = 0.5
     else:
         alpha = alpha
 
-    # If posinf is negative or imaginary set it to its modulus
-    posinf = np.abs(posinf);
+    # If damp is negative or imaginary set it to its modulus
+    damp = np.abs(damp);
 
     # Define initial K-matrices
-    KL = (E + 1j * posinf) * s_l - h_l;
-    K0 = (E + 1j * posinf) * s_0 - h_0;
-    KR = (E + 1j * posinf) * s_r - h_r;
+    KL = (E + 1j * damp) * s_l - h_l;
+    K0 = (E + 1j * damp) * s_0 - h_0;
+    KR = (E + 1j * damp) * s_r - h_r;
 
     # Form initial guess
     AN = K0 - initialguess
@@ -374,3 +372,99 @@ def simple_iterative_greens_function(E, h_l, h_0, h_r, **kwargs):
 
     return NewSelfEnergy
 
+
+def sancho_rubio_iterative_greens_function(E, h_l, h_0, h_r, **kwargs):
+    """Computes the self-energy matrix using the Sancho-Rubio technique.
+    Parameters
+    ----------
+    E : float
+        energy
+
+    h_l : numpy array
+          left block of three-block-diagonal Hamiltonian
+
+    h_0 : numpy array
+          central block of three-block-diagonal Hamiltonian
+
+    h_r : numpy array
+          right block of three-block-diagonal Hamiltonian
+
+    kwargs : dict
+          additional named arguments:
+            s_l : numpy array (Default value = zero matrix)
+                  left block of three-block-diagonal overlap matrix
+
+            s_0 : numpy array (Default value = identity matrix)
+                  central block of three-block-diagonal overlap matrix
+
+            s_r : numpy array (Default value = zero matrix)
+                  right block of three-block-diagonal overlap matrix
+                  
+           damp : float (Default value = 1e-6)
+                  positive infinitesimal a.k.a broadening parameter
+                  that shifts solution off real axis
+
+          nconv : float (Default value = 1e-10)
+                  convergence requirement for iteration using the 2-norm,
+                  also known as the largest singular value.
+
+    Returns
+    -------
+    numpy.matrix
+
+        #Guide for usage.
+
+        A user provides the energy; onsite-, hopping-, and overlap- 
+        matrices; the broadening parameter; and convergence criterion
+        in order to compute the retarded self-energy.
+
+    """
+
+    s_l = kwargs.get('s_l', np.zeros(h_0.shape[0]))
+    s_0 = kwargs.get('s_0', np.identity(h_0.shape[0]))
+    s_r = kwargs.get('s_r', np.zeros(h_0.shape[0]))
+    damp = kwargs.get('damp', 1e-6)
+    nconv = kwargs.get('nconv', 1e-10)
+
+    # If damp is negative or imaginary set it to its modulus
+    damp = np.abs(damp);
+
+    # Define initial K-matrices
+    KL = (E + 1j * damp) * s_l - h_l;
+    K0 = (E + 1j * damp) * s_0 - h_0;
+    KR = (E + 1j * damp) * s_r - h_r;
+
+    #Seed values
+    alpha0 = KR;
+    beta0 = KL;
+    esOld = K0;
+    eOld = esOld;
+    alpha = alpha0; 
+    beta  = beta0;  
+
+    itr = 0;       #initialise iteration check
+    normcheck = 1; #initialise normalisation check 
+    maxiter=100
+
+    while (itr<maxiter)&(normcheck>nconv):
+        itr = itr + 1; 
+
+        lu, piv = linalg.lu_factor(eOld)
+        a = linalg.lu_solve((lu, piv), alpha)
+        b = linalg.lu_solve((lu, piv), beta)
+
+        alphab = alpha.dot(b);    
+        esNew = esOld - alphab;
+        eNew = eOld - beta.dot(a) - alphab;     
+        alpha = alpha.dot(a);
+        beta = beta.dot(b);
+        normcheck = linalg.norm(esNew-esOld, 2)
+        eOld  = eNew ;
+        esOld = esNew;
+        #print(normcheck) #for debugging
+
+
+    lu, piv = linalg.lu_factor(esNew)
+    SelfEnergy = alpha0.dot( linalg.lu_solve( (lu, piv), beta0) )
+
+    return SelfEnergy

test/test_greens_functions.py

@@ -855,10 +855,10 @@ def test_simple_iterative_greens_function():
                                            nconv=nconv)
 
     # The output should then be:
-    value = np.array([[ 0.24781094-0.80987031j, -0.49651996-0.15338929j],
-                      [-0.49651996-0.15338929j,  0.24781094-0.80987031j]])
-
-    np.testing.assert_allclose(ans, value)
+    value = np.array([[ 0.25-0.814841j, -0.5 -0.153404j],
+                      [-0.5 -0.153404j,  0.25-0.814841j]])
+    
+    np.testing.assert_allclose(ans, value, rtol=1e-05)
 
 
 if __name__ == '__main__':