anchor.diff:  do i = 1 , nx
   do j = 1 , ny
       k = gindex(i,j)
  -    write(433,"(2i,10e20.6)"),i,j,sum( abs(waveL(:,k))**2 ),  sum( abs(waveR(:,k))**2 ), &
  +    write(433,"(2i,100e20.6)"),i,j,sum( abs(waveL(:,k))**2 ),  sum( abs(waveR(:,k))**2 ), &
                                     sum( abs(waveL(:,k))**2 ) + sum( abs(waveR(:,k))**2 ),abs(waveL(:,k))**2
   enddo
       write(433,*),""
  
positive.code:  doubleprecision            :: zeros_vector(200)
   doubleprecision            :: a_dx,a_Emin,a_Emax,a_Bz
   integer                    :: no_expected_states
   integer ,parameter         :: nx = 100
   integer ,parameter         :: ny = 15
   doubleprecision,parameter  :: dx = 5.0 ! [nm]
   integer , dimension(nx,ny) :: gindex ! converts local index (i,j) to global index
   integer :: i,j,k,p,ytip
   doubleprecision            :: lead_translation_vec(3),Ef
   doubleprecision            :: Vpot(nx*ny),Vtip
   complex*16,allocatable     :: H0(:,:),G0(:,:),waveR(:,:),waveL(:,:),waveN(:,:)
   complex*16,allocatable     :: invGamma(:,:),dSkm(:,:),tmpSkm(:,:)
   complex*16 :: L(ny,ny),QdagUdag(ny,ny),tmpGU(ny,ny),tmp2GU(ny,ny)
   type(qsmatrix),dimension(2,2) :: smatrix0
   ! Use atomic units in effective band model -> see modunit.f90 for more details
   call modunits_set_GaAs_params()
   a_dx = dx * L2LA ! convert it to atomic units
   a_Bz = BtoAtomicB(0.001D0) ! 1 Tesla to atomic units
  
   ! Initalize system
   call qt%init_system()
  
   ! ----------------------------------------------------------
   ! 1. Create mesh - loop over width and height of the lattice
   ! ----------------------------------------------------------
  
negative.code: 
          total_Tn = total_Tn  + sum(abs(qt%leads(2)%Tnm(:))**2)
      enddo
  !    print"(100e20.6)",ytip*dx,total_Tn
      write(112,"(100e20.6)"),xtip*dx,ytip*dx,total_Tn
  
  
  
  