[ase2-svncheckins] [svn commit /var/www/svn/CamposASE2] r1085 - trunk/ASE/Transport

svn at fysik.dtu.dk svn at fysik.dtu.dk
Tue Oct 28 15:57:14 CET 2008


Author: carstenr
Date: 2008-10-28 15:57:14 +0100 (Tue, 28 Oct 2008)
New Revision: 1085

Added:
   trunk/ASE/Transport/MySelfEnergies.py
Log:
Added an improved version of the IntSelfEnergy module


Added: trunk/ASE/Transport/MySelfEnergies.py
===================================================================
--- trunk/ASE/Transport/MySelfEnergies.py	                        (rev 0)
+++ trunk/ASE/Transport/MySelfEnergies.py	2008-10-28 14:57:14 UTC (rev 1085)
@@ -0,0 +1,480 @@
+import Numeric as num
+
+def GetHartreeAll(D, V_ijkl):
+    if type(D) == list:
+        D = D[0] + D[1]
+    else:
+        D = 2 * D
+
+    N = len(D)
+    V_H = num.empty([N, N], num.Complex)
+    for i in range(N):
+        for j in range(N):
+            V_H[i, j] = num.dot(num.ravel(V_ijkl[i, :, j, :]), D.flat)
+    return V_H
+
+def GetHartreePartial(D, V_ijij, V_ijji, V_iijj, V_iiij, V_ikjk=None):
+    if type(D) == list:
+        D = D[0] + D[1]
+    else:
+        D = 2 * D
+
+    N = len(D)
+    V_H = num.empty([N, N], num.Complex)
+    if V_iijj is None:
+        V_iijj = num.zeros([N, N], num.Complex)
+    if V_iiij is None:
+        V_iiij = num.zeros([N, N], num.Complex)
+    for i in range(N):
+        for j in range(N):
+            if i == j:
+                if V_ikjk is not None:
+                    V_H[i, i] = num.sum((D * V_ikjk[:, :, i]).flat)
+                else:
+                    V_H[i, i] = (num.dot(num.diagonal(D), V_ijij[i]) +
+                                 num.dot(D[i], V_iiij[i]) +
+                                 num.dot(D[:, i], num.conjugate(V_iiij[i])) -
+                                 2 * D[i, i] * V_iiij[i, i])
+            else:
+                V_H[i, j] = (D[i, j] * V_iijj[i, j] +
+                             D[j, i] * V_ijji[i, j] +
+                             D[i, i] * V_iiij[i, j] +
+                             D[j, j] * num.conjugate(V_iiij[j, i]))
+                if V_ikjk is not None:
+                    V_H[i, j] += (num.dot(num.diagonal(D), V_ikjk[i, j, :]) -
+                                  D[i, i] * V_iiij[i, j] -
+                                  D[j, j] * num.conjugate(V_iiij[j, i]))
+    return V_H
+
+def GetFockAll(D, V_ijkl):
+    if type(D) == list:
+        return [GetFockAll(D[0], V_ijkl), GetFockAll(D[1], V_ijkl)]
+
+    N = len(D)
+    V_F = num.empty([N, N], num.Complex)
+    for i in range(N):
+        for j in range(N):
+            V_F[i, j] = -num.dot(num.ravel(V_ijkl[i, :, :, j]), D.flat)
+    return V_F
+
+def GetFockPartial(D, V_ijij, V_ijji, V_iijj, V_iiij, V_ikjk=None):
+    if type(D) == list:
+        return [GetFockAll(D[0], V_ijkl), GetFockAll(D[1], V_ijkl)]
+    
+    N = len(D)
+    V_F = num.empty([N, N], num.Complex)
+    if V_iijj is None:
+        V_iijj = num.zeros([N, N], num.Complex)
+    if V_iiij is None:
+        V_iiij = num.zeros([N, N], num.Complex)
+    for i in range(N):
+        for j in range(N):
+            if i == j:
+                V_F[i, i] = -(num.dot(num.diagonal(D), V_ijji[i]) +
+                              num.dot(D[i], V_iiij[i]) +
+                              num.dot(D[:, i], num.conjugate(V_iiij[i])) -
+                              2 * D[i, i] * V_iiij[i, i] +
+                              D[i, i] * V_ijij[i, i] - D[i, i] * V_ijji[i, i])
+            else:
+                V_F[i, j] = -(D[j, i] * V_ijij[i, j] +
+                              D[i, j] * V_iijj[i, j] +
+                              D[i, i] * V_iiij[i, j] +
+                              D[j, j] * num.conjugate(V_iiij[j, i]))
+                if V_ikjk is not None:
+                    V_F[i, j] -= (num.dot(D[:, i], V_ikjk[:, j, i]) -
+                                  D[i,i]*V_ikjk[i,j,i] - D[j,i]*V_ikjk[j,j,i] +
+                                  num.dot(D[j, :], V_ikjk[i,:,j]) -
+                                  D[j,j]*V_ikjk[i,j,j] - D[j,i]*V_ikjk[i,i,j])
+    return V_F
+
+class NonEqIntSelfEnergy:
+    def __init__(self, nonintgreenfunction, name):
+	self.SetNonIntGreenFunction(nonintgreenfunction)
+	self.SetActive(True)
+	self.SetName(name)
+	self.shape = nonintgreenfunction.GetMatrixDimension()
+
+    def SetName(self,name):
+	self.name=name
+
+    def GetName(self):
+	return self.name
+
+    def Correlated(self):
+	pass
+
+    def SetNonIntGreenFunction(self, nonintgf):
+	self.nonintgf = nonintgf
+
+    def GetNonIntGreenFunction(self):
+	return self.nonintgf
+
+    def GetRetarded(self, index, conv=False):
+	# Always use this method for the retarded self energy (don't ask
+	# for the list).
+	# This ensures proper updates. 
+	# Remember to ask for Update() before returning the self-energy. 
+	pass
+
+    def GetLesser(self, index):
+	# Always use this method for the retarded self energy (don't ask
+	# for the list).
+	# This ensures proper updates. 
+	# Remember to ask for Update() before returning the self-energy. 
+	pass
+
+    def Active(self):
+        return self.activity
+
+    def SetActive(self, activity):
+	self.activity = activity
+
+    def Unregister(self):
+	pass
+
+    def Update(self, conv=False):
+	pass
+        
+class NonEqConstantSelfEnergy(NonEqIntSelfEnergy):
+    def __init__(self, nonintgreenfunction, sigma=None, name='constant'):
+	NonEqIntSelfEnergy.__init__(self, nonintgreenfunction, name)
+	self.SetSigma(sigma)
+
+    def GetSigma(self):
+	if not hasattr(self, 'sigma'):
+	    self.Update()
+	return self.sigma
+
+    def SetSigma(self, sigma):
+	self.sigma = sigma
+
+    def Correlated(self):
+	return False
+
+    def GetLesser(self,index):
+	return num.zeros(self.shape, num.Complex)
+
+    def GetGreater(self,index):
+	return num.zeros(self.shape, num.Complex)
+	
+    def GetRetarded(self,index,converged=False):
+	return self.GetSigma()   
+
+    def Unregister(self):
+	if hasattr(self, 'sigma'):
+	    del self.sigma
+
+class NonEqHartreeSelfEnergy(NonEqConstantSelfEnergy):
+    def __init__(self, nonintgreenfunction, Vlst, initialhartree=None):
+	NonEqConstantSelfEnergy.__init__(self, nonintgreenfunction,
+					 name='hartree')
+        N = self.shape[0]
+	while len(Vlst) < 4:
+	    Vlst.append(num.zeros((N, N), num.Complex))
+        if len(Vlst) < 5:
+            Vlst.append(None)
+	self.vlst = Vlst
+
+	# Set the hartree potential already included in H0
+	self.initialhartree = initialhartree
+	self.Update()
+
+    def GetInitialHartree(self):
+	if self.initialhartree == None:
+	    print "###########################################################"
+	    print "### OBS: Setting initial Hartree potential to           ###"
+            print "###      equilibrium Hartree potential                  ###"
+	    print "###########################################################"
+            mul = self.GetNonIntGreenFunction().GetLeftChemicalPotential()
+            mur = self.GetNonIntGreenFunction().GetRightChemicalPotential()
+            assert abs(mul - mur) < 0.00001, "ERROR: You must set the chemical potentials to zero before initializing Hartree potential"
+	    self.initialhartree = self.GetFullHartree()
+	return self.initialhartree
+
+    def GetVlst(self):
+	return self.vlst
+
+    def GetFullHartree(self):
+	v = self.GetVlst()
+	D = self.GetNonIntGreenFunction().GetDensityMatrix()
+	return GetHartreePartial(D, v[0], v[1], v[2], v[3], v[4])
+
+    def Update(self, conv=False):
+	self.SetSigma(self.GetFullHartree()-self.GetInitialHartree())
+
+class NonEqFockSelfEnergy(NonEqConstantSelfEnergy):
+    def __init__(self, nonintgreenfunction, Vlst, Fcore=None):
+	NonEqIntSelfEnergy.__init__(self, nonintgreenfunction, 'fock')
+        N = self.shape[0]
+	while len(Vlst) < 4:
+	    Vlst.append(num.zeros((N, N), num.Complex))
+        if len(Vlst) < 5:
+            Vlst.append(None)
+	self.vlst = Vlst
+        self.Fcore = Fcore
+
+    def GetVlst(self):
+	return self.vlst
+
+    def SetSigma(sigma):
+        if self.Fcore == None:
+            return NonEqConstantSelfEnergy.SetSigma(self, sigma)
+        return NonEqConstantSelfEnergy.SetSigma(self, sigma + self.Fcore)
+
+    def Update(self,conv=False):
+	v = self.GetVlst()
+	D = self.GetNonIntGreenFunction().GetDensityMatrix()
+        self.SetSigma(GetFockPartial(D, v[0], v[1], v[2], v[3], v[4]))
+        if self.Fcore != None:
+            self.sigma += Fcore
+
+class HartreeAll(NonEqHartreeSelfEnergy):
+    """This takes a V list, with the single element V_ijkl"""
+    def GetFullHartree(self):
+	V = self.GetVlst()[0]
+	D = self.GetNonIntGreenFunction().GetDensityMatrix()
+	return GetHartreeAll(D, V)
+
+class FockAll(NonEqFockSelfEnergy):
+    """This takes a V list, with the single element V_ijkl"""
+    def Update(self, conv=False):
+	V = self.GetVlst()[0]
+        D = self.GetNonIntGreenFunction().GetDensityMatrix()
+        self.SetSigma(GetFockAll(D, V))
+    
+class HartreeAtomAll(NonEqHartreeSelfEnergy):
+    """This takes a V list, with the elements V_ijij, V_ijji, V_iijj, V_iiij
+    and 'basis', which is of the form
+      [(index of first basis func for first atom, last basis for first atom),
+       (first basis of 2. atom, last basis of 2. atom),
+       ...]
+    and Va, which is a list of the V_ijkl matrix for each atom
+    """
+    def __init__(self, g0, Vlst, basis, Va, initialhartree=None):
+        self.basis = basis
+        self.Va = Va
+        NonEqHartreeSelfEnergy.__init__(self, g0, Vlst, initialhartree)
+        
+    def GetFullHartree(self):
+        V_H = NonEqHartreeSelfEnergy.GetFullHartree(self)
+        for (start, end), V in zip(self.basis, self.Va):
+            D = self.GetNonIntGreenFunction().GetDensityMatrix()[
+                start:end, start:end].copy()
+
+            v1 = self.vlst[0][start:end, start:end]
+            v2 = self.vlst[1][start:end, start:end]
+            v3 = self.vlst[2][start:end, start:end]
+            v4 = self.vlst[3][start:end, start:end]
+            v5 = None
+            if self.vlst[4] is not None:
+                v5 = self.vlst[4][start:end, start:end, start:end]
+            
+            V_H[start:end, start:end] += GetHartreeAll(D, V)
+            V_H[start:end, start:end] -= GetHartreePartial(D,v1,v2, v3, v4, v5)
+        return V_H
+
+class FockAtomAll(NonEqFockSelfEnergy):
+    """This takes a V list, with the elements V_ijij, V_ijji, V_iijj, V_iiij
+    and 'basis', which is of the form
+      [(index of first basis func for first atom, last basis for first atom),
+       (first basis of 2. atom, last basis of 2. atom),
+       ...]
+    and Va, which is a list of the V_ijkl matrix for each atom
+    """
+    def __init__(self, g0, Vlst, basis, Va):
+        self.basis = basis
+        self.Va = Va
+        NonEqFockSelfEnergy.__init__(self, g0, Vlst)
+        
+    def Update(self, conv=False):
+        NonEqFockSelfEnergy.Update(self, conv=False)
+        V_F = self.GetSigma()
+        for (start, end), V in zip(self.basis, self.Va):
+            D = self.GetNonIntGreenFunction().GetDensityMatrix()[
+                start:end, start:end].copy()
+
+            v1 = self.vlst[0][start:end, start:end]
+            v2 = self.vlst[1][start:end, start:end]
+            v3 = self.vlst[2][start:end, start:end]
+            v4 = self.vlst[3][start:end, start:end]
+            v5 = None
+            if self.vlst[4] is not None:
+                v5 = self.vlst[4][start:end, start:end, start:end]
+            
+            V_F[start:end, start:end] += GetFockAll(D, V)
+            V_F[start:end, start:end] -= GetFockPartial(D, v1, v2, v3, v4, v5)
+        self.SetSigma(V_F)
+
+class PhononHartreeSelfEnergy(NonEqConstantSelfEnergy):
+    def __init__(self, nonintgf, M_lii, Omega_l):
+        NonEqConstantSelfEnergy.__init__(self, nonintgf, name='phonon_hartree')
+
+        # Set the parameters. This also forces an Unregister
+        self.SetMassMatrix(M_lii)
+        self.SetFrequencyList(Omega_l)
+
+    def GetSigma(self):
+        return self.GetRetarded(0)
+
+    def GetRetarded(self, energy, converged=False):
+        if self.update_needed:
+            self.Update()
+        return self.retarded_ii
+
+    def SetMassMatrix(self, M_lii):
+        self.M_lii = M_lii
+        self.Unregister()
+        
+    def SetFrequencyList(self, Omega_l):
+        self.Omega_l = Omega_l
+        self.Unregister()
+
+    def Update(self, conv=False):
+        """
+          -- H,r  --     2 i     l  /    --  <     l
+          >     = >   --------- M   |dw  >  G (w) M
+          --      --  pi Omega   ij /    --  mn  nm
+            ij    l           l          mn
+        """
+        # Allocate space for retarded array
+        self.retarded_ii = num.zeros(self.shape, num.Complex)
+
+        # Get lesser Green function
+        less_eii = self.nonintgf.GetLesserList()
+
+        # Energy grid-spacing
+        de = self.nonintgf.GetEnergyList()[1] -self.nonintgf.GetEnergyList()[0]
+
+        # Loop over eigen-frequencies
+        for Omega, M_ii in zip(self.Omega_l, self.M_lii):
+            # Integrate over the trace of G_less times M
+            ## int_trace = num.sum(num.trace(num.dot(
+            ##     less_eii, M_ii), axis1=1, axis2=2)) * de
+            int_trace = 0.
+            for less_ii in less_eii: # loop necesary because of num.dot bug
+                int_trace += num.trace(num.dot(less_ii, M_ii))
+            int_trace *= de
+            
+            # Update retarded phonon Hartree self-energy
+            self.retarded_ii += M_ii * 2.j * int_trace / (num.pi * Omega)
+
+        self.update_needed = False
+
+    def Unregister(self):
+        self.update_needed = True
+
+class PhononFockSelfEnergy(NonEqIntSelfEnergy):
+    def __init__(self, nonintgf, M_lii, Omega_l):
+        NonEqIntSelfEnergy.__init__(self, nonintgf, 'phonon_fock')
+
+        # Set the parameters. This also forces an Unregister
+        self.SetMassMatrix(M_lii)
+        self.SetFrequencyList(Omega_l)
+
+    def GetRetarded(self, index, converged=False):
+        if self.update_needed:
+            self.Update()
+        return self.retarded_eii[index]
+
+    def GetGreater(self, index):
+        if self.update_needed:
+            self.Update()
+        return self.greater_eii[index]
+
+    def GetLesser(self, index):
+        if self.update_needed:
+            self.Update()
+        return self.lesser_eii[index]
+
+    def SetMassMatrix(self, M_lii):
+        self.M_lii = M_lii
+        self.Unregister()
+        
+    def SetFrequencyList(self, Omega_l):
+        self.Omega_l = Omega_l
+        self.Unregister()
+
+    def Update(self, conv=False):
+        """Update the retarded, greater, and lesser self-energies according to
+        the rules below, where it has been assumed that the boson temperature
+        is zero::
+        
+          -- f,r      -- --  l  /     l               l    \  l
+          >     (w) = >  >  M   | diff (w) - i H[diff] (w) | M
+          --          -- --  im \     mn              mn   /  nj
+            ij        l  mn       
+
+                    l      1 /  >                <             \
+          where diff (w) = - | G (w - Omega ) - G (w + Omega ) |
+                    ij     2 \  ij         l     ij         l  /
+
+          and H[diff] is the Hilbert transform of this (with the sign
+          convention H[f](y) = 1/pi p.v. int dx f(x) / (x - y) ).
+
+          -- f,<        -- -- l  <              l
+          >     (w)  =  >  > M  G (w + Omega ) M
+          --            -- -- im mn         l   nj
+            ij          l  mn    
+
+          -- f,>        -- -- l  >              l
+          >     (w)  =  >  > M  G (w - Omega ) M
+          --            -- -- im mn         l   nj
+            ij          l  mn    
+                        
+        """
+        from ASE.Utilities.ArrayTools import TranslateAlongAxis0NonPeriodic as Translate
+
+        # Import tools for making Hilbert transforms
+        from ASE.Transport.Hilbert import Hilbert
+
+        # Allocate space for retarded, greater, and lesser arrays
+        energies_e = self.nonintgf.GetEnergyList()
+        Ne = len(energies_e)
+        de = energies_e[1] - energies_e[0]
+        dim = (Ne,) + self.shape
+        self.retarded_eii = num.zeros(dim, num.Complex)
+        self.greater_eii  = num.zeros(dim, num.Complex)
+        self.lesser_eii   = num.zeros(dim, num.Complex)
+
+        # Get lesser and greater Green functions
+        Ggreat_eii = self.nonintgf.GetGreaterList()
+        Gless_eii  = self.nonintgf.GetLesserList()
+
+        # Temporary arrays
+        min_eii  = num.zeros(dim, num.Complex)
+        plus_eii = num.zeros(dim, num.Complex)
+
+        # Loop over eigen-frequencies
+        for Omega, M_ii in zip(self.Omega_l, self.M_lii):
+            # Translation steps
+            steps = Omega / de
+
+            # Interpolate if steps is non-integer
+            dec = steps % 1
+            steps = int(steps - dec)
+            min_eii.flat[:] = 0.
+            plus_eii.flat[:] = 0.
+            for weight in (1 - dec, dec):
+                if round(weight, 5) != 0:
+                    min_eii += weight * Translate(Ggreat_eii, +steps,
+                                                  num.Complex)
+                    plus_eii += weight * Translate(Gless_eii, -steps,
+                                                      num.Complex)
+                steps += 1
+
+            # Make Hilbert transform
+            hilbert_eii = Hilbert(min_eii - plus_eii)
+            
+            # Update retarded, greater, and lesser self-energies at each energy
+            for e in range(Ne): # loop necesary because of num.dot bug
+                self.retarded_eii[e] += num.dot(M_ii, num.dot(
+                    .5 * (min_eii[e] - plus_eii[e]) - .5j * hilbert_eii[e]
+                    , M_ii))
+
+                self.greater_eii[e] += num.dot(M_ii, num.dot(min_eii[e], M_ii))
+                self.lesser_eii[e] += num.dot(M_ii, num.dot(plus_eii[e], M_ii))
+        
+        self.update_needed = False
+
+    def Unregister(self):
+        self.update_needed = True



More information about the ase2-svncheckins mailing list