subroutine computeSelfEnergy(energy,N_a,N_l_axialcell,HoppingMatrixArray,LayerHamiltonian, &
& eta,routineSelector,N_iterations,SL,SR)
! by Luis E. F. Foa-Torres (luisfoa@gmail.com), October 2005. 
! This example code is provided 'as it is' and only for educational purposes.

! An ordered system is assumed in the calculations.

! N_iterations: number of iterations in the decimations procedure - equivalent to 
!the decimation of 2^(N_iterations-1) unit cells in the axis direction.

! HoppingMatrixArray(1:N_l_axialcell,1:N_a,1:N_a): array containing the hopping matrices for a unit cell in the axial direction.
!
!For Armchair Graphene Nanoribbons, N_l_axialcell is usually set to 4 (i.e. a unit cell of 4 layers), 
!and HoppingMatrixArray contains the matrices, beta_1, beta_2, 
!transpose(beta_1), transpose(beta_2), this is, the interlayer hoppings within the unit cell.
!
!On output, the routine gives the self-energy correction to the left (SL) or to the right (SR) of a given layer.
!

!N_l_axialcell: number of layers in the axial direction contained in a "unit cell".
!N_a: number of atoms in each layer (circumferential direction).
!eta: regularization parameter, usually small, should be tuned for the value of N_iterations in use.


  implicit none
  complex*8, parameter  :: cIm=(0.d0,1.d0)
  integer, intent(in)   :: N_a, N_l_axialcell, N_iterations, routineSelector
  integer               :: temp_index, condnumber
  integer               :: i, ii, jj      !loop variables.
  real*8, intent(in)    :: energy
  complex*16, intent(in) :: eta
  complex*16, intent(in) :: HoppingMatrixArray(1:N_l_axialcell,1:N_a,1:N_a)
  complex*16, intent(in) :: LayerHamiltonian(1:N_a,1:N_a)
  complex*16, intent(out):: SL(1:N_a,1:N_a), SR(1:N_a,1:N_a)
  complex*16 :: IdMatrix(1:N_a,1:N_a) 
!  complex*16 :: V_L(1:N_a,1:N_a), V_R(1:N_a,1:N_a)
  complex*16 :: LayerG(1:N_a,1:N_a), contribSL(1:N_a,1:N_a), contribSR(1:N_a,1:N_a)
  complex*16 :: Veff(1:N_a,1:N_a), VL(1:N_a,1:N_a), VR(1:N_a,1:N_a)
  complex*16 :: Vefft(1:N_a,1:N_a), VLt(1:N_a,1:N_a), VRt(1:N_a,1:N_a)

  !Initialize decimation variables
  SL=0.d0
  SR=0.d0

  call ComplexIdentityMatrix(N_a,IdMatrix)
  !Define Hoppings (right to left)
  do jj=1,N_a
    do ii=1,N_a
      VL(ii,jj)=HoppingMatrixArray(1,ii,jj)
      temp_index=2-int(1/N_l_axialcell)*N_l_axialcell  !CHECK!!!!
      VR(ii,jj)=HoppingMatrixArray(temp_index,ii,jj)
    end do
  end do
!!!!
  VLt=transpose(VL)
  VLt=conjg(VLt)
  VRt=transpose(VR)
  VRt=conjg(VRt)
!!!!
  !Decimate internal layers from left to right.

  do i=1, N_l_axialcell-1
    LayerG=(energy+cIm*eta)*IdMatrix-LayerHamiltonian-SR
    call Invertila(LayerG, N_a, routineSelector, condnumber)     
    contribSL=0.d0; contribSR=0.d0; Veff=0.d0
    call decimateLayer(N_a, LayerG, VL, VLt, VR, VRt, contribSL, contribSR, Veff, Vefft)
    !  Sigma_L=V_L*LayerG*conjg(V_L)
    !  Sigma_R=V_R*LayerG*conjg(V_R)
    !  Veff=V_L*LayerG*V_R
   
    SL=SL+contribSL
    SR=contribSR    
    VL=Veff
    VLt=Vefft
    !!!!Update V_R for the next iteration in the decimation procedure.!!!!
    temp_index=i+2-int((i+1)/N_l_axialcell)*N_l_axialcell  
    !This gives the series: 1,2,... N_l_axialcell, 1,2...
    do jj=1,N_a
      do ii=1,N_a
        VR(ii,jj)=HoppingMatrixArray(temp_index, ii,jj)
      end do
    end do
    VRt=transpose(VR)
    VRt=conjg(VRt) 
  end do

!Once the internal layers in the (axial) unit cell have been decimated we have an 
!effective system where all the sites are equivalent.
!Now we start decimating the even layers there until convergence of the 
!self-energy corrections is achieved. 
!Note that instead of decimating the layers one by one we do that with all the 
!even layers at the same time. Since the computational cost is the same in both 
!cases, we boost the convergence speed.

  !Define the LayerHamiltonian and InterLayerHoppings.
  !Decimate even layers from left to right. After N steps we have decimated a 
  !system of length 2^(N-1).
  if(N_l_axialcell .ne. 1) then
    VL=Veff; VR=Veff
    VLt=Veff; VRt=Veff 
  end if

  
  do i=1, N_iterations
    LayerG=(energy+cIm*eta)*IdMatrix-LayerHamiltonian-(SR+SL)
    call Invertila(LayerG, N_a, routineSelector, condnumber)     
    contribSL=0.d0; contribSR=0.d0; Veff=0.d0
    
    call decimateLayer(N_a, LayerG, VL, VLt, VR, VRt, contribSL, contribSR, Veff, Vefft)
   
    SL=SL+contribSL
    SR=SR+contribSR     ! Note that in this case the expressions are symmetric 
    !because we are decimating sites everywhere and not only from left to right.
    VL=Veff; VLt=Vefft
    !!!!Update V_R for the next iteration in the decimation procedure.!!!!
    VR=Veff; VRt=Vefft
!    VR=conjg(VR)
  end do

  ! Output variables: Sigma_L=SL; Sigma_R=SR

end subroutine computeSelfEnergy



subroutine Invertila(A, N, routineSelector, condnumber)
! Computes the inverse of the square matrix A(1:N, 1:N)
! The routine used to this end depends on the value of the integer routineSelector.
! Uses Lapack

    INTEGER            INFO, LDA, LWORK, N, routineSelector
    INTEGER            IPIV(1:N), condnumber
    COMPLEX*16      :: A(1:N,1:N), Work(1:N) ! double complex

	if(routineSelector .eq. 1) then  ! Use Lapack, add appropriate options when compiling
      LDA=N
      LWork=N
      CALL ZGETRF( N, N, A, LDA, IPIV, INFO)
      CALL ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) !Lapack
!    else ! Else you can add your own routine below and uncomment this line.
!      CALL cInv(A, N, condnumber)! G = (E_Heff)^(-1) 
    end if    
end subroutine Invertila
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Subroutine defineLayerHamiltonian defines the Layer and InterLayer matrices.
! This sample routine works for ArmchairGrapheneNanoribbons only

! On output gives the matrices LayerBlock and InterLayerHopping.
!#################################
!Note on ussage:
!Note that the interlayerHopping is not the same from one layer to the next one, but every two layers.
!The interLayerHopping follows the sequence: beta_1,beta_2, Transpose_beta_1, Transpose_beta_2,...
!#################################
Subroutine defineLayerHamiltonian_ACGNR(N_a, E_o, V, LayerBlock, beta_1, beta_2)
!N_a : number of A atoms in each layer.
!For the present case we assume that N_b=N_a, so that N_c=total number of atoms in the selected layer.
!U_a (U_b): potential energy at the A (B) atoms.
!V : C-C hopping matrix element (pi orbitals).
  integer :: N_a
  integer :: j
  real*8      :: Id_N_a(1:N_a,1:N_a)
  complex*16  :: E_o, V
  complex*16  :: beta_1(1:N_a,1:N_a), beta_2(1:N_a,1:N_a), tempArray_N_a(1:N_a,1:N_a)
  complex*16  :: LayerBlock(1:N_a,1:N_a)
  

  !define the matrices corresponding to the intra and inter-layer coupling. 
  !
  !definitions valid for ARMCHAIR GNRs.
  !

  call IdentityMatrix(N_a,Id_N_a)
  LayerBlock=E_o*Id_N_a
  beta_2= V * Id_N_a
  beta_1= V * Id_N_a
  
  do j=1,N_a-1
    beta_1(j,j+1)=V
  end do
!  beta_2B(N_a,1)=V
 
end subroutine defineLayerHamiltonian_ACGNR