Neq

hubbard/hamiltonian.py

@@ -65,9 +65,12 @@ def __init__(self, TBHam, n=0, U=None, Uij=None, q=(0., 0.), nkpt=[1, 1, 1], kT=
         H0 = self.TBHam.copy()
         H0.shift(np.pi) # Apply a shift to incorporate effect of S
         if self.spin_size > 1:
-            s = H0.H.tocsr(0).data.tostring() + H0.H.tocsr(1).data.tostring()
+            #s = H0.H.tocsr(0).data.tostring() + H0.H.tocsr(1).data.tostring()
+            s = H0.H.tocsr(0).data + H0.H.tocsr(1).data
+
         else:
-            s = H0.H.tocsr(0).data.tostring()
+            #s = H0.H.tocsr(0).data.tostring()
+            s = H0.H.tocsr(0).data
         self._hash_base = s
         del H0
 

hubbard/negf.py

@@ -211,7 +211,7 @@ def convert2SelfEnergy(elec_SE, mu):
                             # And for each point in the Neq CC
                             self._cc_neq_SE[spin][ik][ic][i] = se.self_energy(cc, k=k, **kw)
 
-    def calc_n_open(self, H, q, qtol=1e-5, a=1.):
+    def calc_n_open(self, H, q, qtol=1e-5, method='residual', gamma=1.0,lambda_tf=7.4, self_capacitance=True, C_self = 0.1, mixing=None,alpha=0.2):
         """
         Method to compute the spin densities from the non-equilibrium Green's function
 
@@ -224,6 +224,15 @@ def calc_n_open(self, H, q, qtol=1e-5, a=1.):
         qtol: float, optional
             tolerance to which the charge is going to be converged in the internal loop
             that finds the potential of the device (i.e. that makes the device neutrally charged)
+        method: str
+            model in which the charge will be calculated. The options are residual 'residual' or thomas fermi 'tf'
+        lambda_tf: float
+            Thomas Fermi wavelength from https://arxiv.org/pdf/1010.0508 1.42 Å
+            Appl. Phys. Lett. 92, 123110 (2008)                          7.4  Å
+        C_self: str
+            Self-capacitance term, e per eV per atom (i.e., Kii = 10 eV/e)
+        mixing: str, optional
+            	Use 'damped' for improve convergence with dq
 
         Returns
         -------
@@ -363,8 +372,7 @@ def calc_n_open(self, H, q, qtol=1e-5, a=1.):
 
         # Save Fermi-level of the device
         self.Ef = Ef
-
-
+
         if self.NEQ:
             # Add potential in each site depending on how much the neq charges deviate from the eq situation. This term comes from the extended Huckel model
             # a should be positive: if q_neq > q_eq then the potential should rise (less favourable for electrons to occupy that site)
@@ -373,7 +381,32 @@ def calc_n_open(self, H, q, qtol=1e-5, a=1.):
             q_eq = self.H_eq.n.sum(axis=0)
             dq = (q_neq - q_eq)[self.a_dev]
             E = self.H.TBHam.tocsr(spin).diagonal()[self.a_dev]
-            self.H.TBHam[self.a_dev, self.a_dev] = E + a * dq
+
+            # Modification by Alan Anaya
+
+            if method == 'residual':
+                self.H.TBHam[self.a_dev, self.a_dev] = E + gamma * dq
+
+            if method == 'tf':
+                # Loop over atoms in the Device
+            	for i,atoms_i in enumerate(self.a_dev):
+            	    #  On-site Change on Atom i
+            	    dV = 0
+            	    for j ,atoms_j in enumerate(self.a_dev):
+            	    	if atoms_i != atoms_j:
+                	        rij	= self.H.TBHam.geometry.rij(atoms_i,atoms_j)
+                	        Kij = (14.4/rij)*np.exp(-rij/lambda_tf) #  1/(4πε0) in eV/Å
+                	        dV += Kij*dq[j]
+            	    # Self-capacitance term
+            	    if self_capacitance==True:
+                	    Kii = 1.0 / C_self
+                	    dV += Kii * dq[i]
+            	    # For debugging
+                    #print(atoms_i,dV)
+            	    if mixing == None:
+            	        self.H.TBHam[atoms_i,atoms_i] += dV
+            	    if mixing == 'damped':
+            	        self.H.TBHam[atoms_i,atoms_i] += alpha*dV
         # Return spin densities and total energy, if the Hamiltonian is not spin-polarized
         # multiply Etot by 2 for spin degeneracy
         return ni, (2./H.spin_size)*Etot

tests_NEQ/0V_gate.py

@@ -0,0 +1,113 @@
+import sisl
+from hubbard import HubbardHamiltonian, sp2, density, plot, NEGF
+import ase
+from ase.visualize import view
+import tk
+import matplotlib.pyplot as plt
+from ase.build import graphene_nanoribbon
+import numpy as np
+import netCDF4
+from sisl import geom
+import os
+
+# Set Hubbard U 
+U=3.5
+kx=1
+
+newpath = f'./Eq_gate_{U}U_{kx}kx' 
+if not os.path.exists(newpath):
+    os.makedirs(newpath)
+
+# Read Device
+device=sisl.get_sile('device.fdf').read_geometry()
+# Remove H from device
+device=device.remove(device.atoms.Z==1)
+device.set_nsc([3,1,1])
+
+# Read Electrode
+elec=sisl.get_sile('elec.xyz').read_geometry()
+elec.set_nsc([3,1,3])
+elec_ts=elec.tile(2,axis=2)
+
+# get unbias hamiltonian
+H0 = sp2(elec_ts, t1=2.7, t2=0.0, t3=0.0, spin='polarized')
+q0 = H0.q0
+
+for gate in np.round(np.arange(-0.25,0.375,0.125),3):
+    #############################################################Electrodes#######################################################################
+    # calcullate charge on each atom
+	dn=gate/q0
+    #assign charge to electrode hamiltonian
+	H_elec = sp2(elec_ts, t1=2.7, t2=0.0, t3=0.0, spin='polarized',dq=dn)
+    #get charge on the gated electrode
+	q=H_elec.q0
+    #pass equal amount of charge on MFH object
+	MFH_elec = HubbardHamiltonian(H_elec, U=U, nkpt=[30,1,30],q=(0.5*q, 0.5*q), kT=0.025) 
+	MFH_elec.random_density()
+	MFH_elec.converge(density.calc_n, tol=1e-6, print_info=True, steps=10)
+	print('Fermi level of'+str(MFH_elec.fermi_level()))
+	MFH_elec.shift(MFH_elec.fermi_level())
+	    #############################################################Plot Bands#######################################################################
+	fig, ax = plt.subplots(figsize=(2,8))
+	fig.suptitle(f'Spin polarization e {gate}',size=24)
+	bs = sisl.BandStructure(MFH_elec.H, [[0.5, 0, 0],[0, 0, 0], [0.0, 0, 0.5]], 100,['X', r'Gamma', 'Z'])
+	lk = bs.lineark(ticks=False)
+	bs_eig = bs.apply.array.eigh()
+	plt.plot(lk, bs_eig,c='red')
+	plt.ylim([-4,4])
+	plt.fill_between(lk, y1=-10, y2=0,color='red',alpha=0.25)
+	plt.xlabel('k', size=14)
+	plt.ylabel('E (eV)', size=14)
+	fig.savefig(f'{newpath}/Electrode_Bandstructure_{gate}.png',dpi=200,bbox_inches="tight")
+	plt.clf()
+	#############################################################Device#######################################################################
+
+    # calculate charge on device on each atom 
+	dQ=(gate)*8
+	dN=dQ/device.na
+	print(dN,dQ)
+	HC= sp2(device, t1=2.7, t2=0.0, t3=0.0, spin='polarized',dq=dN)
+	# with periodic boundary conditions
+	HC.set_nsc([3,1,1])
+    # get total charge
+	Q=HC.q0
+	# MFH object of the central region, same kT and U as electrodes!
+	MFH_HC = HubbardHamiltonian(HC, U=U, kT=0.025, nkpt=[kx,1,1],q=(0.5*Q, 0.5*Q))
+	MFH_HC.set_polarization([112,120,128,138,147,152], dn=[115,121,129,139,146,155])
+	elec_indx = [range(len(H_elec)), range(-len(H_elec), 0)]
+	negf = NEGF(MFH_HC, [(MFH_elec, '-C'), (MFH_elec, '+C')], elec_indx)
+	MFH_HC.converge(negf.calc_n_open, steps=5, tol=1e-4, print_info=True)
+	    #############################################################Plot SP#######################################################################
+	ds=MFH_HC.n[0]-MFH_HC.n[1]
+	dE=MFH_HC.H.tocsr(dim=0).diagonal()-MFH_HC.H.tocsr(dim=1).diagonal()
+	# Create just a figure and only one subplot
+	fig, ax = plt.subplots(nrows=2,ncols=1,figsize=(10,12),sharex=True)
+	fig.suptitle(f'Spin Polarization e {dQ}',size=24)
+
+	#spin polarization
+	pcm1=ax[0].scatter(device[:,2],device[:,0],c=ds,cmap='RdBu_r',vmin=-0.25,vmax=0.25,edgecolors='k',s=300)
+	fig.colorbar(pcm1, ax=ax[0],).set_label(label=r'${Q_\uparrow}-{Q_\downarrow}$',size=15,weight='bold')
+	ax[0].set_ylabel('x Å',size=16)
+	for i, txt in enumerate(ds):
+		if np.abs(txt)>0.05:
+			ax[0].annotate(np.round(txt,2), (device.xyz[i,2]-1.8, device.xyz[i,0]),size=10)
+	# onsite difference        
+	pcm2 = ax[1].scatter(device[:,2],device[:,0],c=dE,cmap='RdBu_r',vmin=-0.5,vmax=0.5,edgecolors='k',s=300)
+	fig.colorbar(pcm2, ax=ax[1]).set_label(label=r'${E_\uparrow}-{E_\downarrow}$',size=15,weight='bold')
+	ax[1].set_ylabel('x Å',size=16)
+	ax[1].set_xlabel('z Å',size=16)
+	fig.tight_layout()
+	fig.savefig(f'{newpath}/SP_{dQ}.png',dpi=200,bbox_inches="tight")
+	plt.clf()
+
+###########################################################Plot DOS######################################################################
+	fig, ax = plt.subplots(figsize=(2,8),)
+	DOS_0=MFH_HC.DOS(np.arange(-1,1,0.01),spin=0)
+	DOS_1=MFH_HC.DOS(np.arange(-1,1,0.01),spin=1)
+	plt.plot(DOS_0,np.arange(-1,1,0.01),label='spin0',c='r')
+	plt.plot(DOS_1,np.arange(-1,1,0.01),label='spin1',c='b')
+	plt.ylabel('E (eV)',size=16)
+	plt.xlabel('DOS',size=16)
+	fig.savefig(f'{newpath}/DOS_{dQ}.png',dpi=200,bbox_inches="tight")
+	plt.clf()
+	print("finish calculation")