anchor.diff:  call qt%save_system(output_folder//"system.xml")
   
   ! Solve scattering problem for Ef=0.001
  -Ef = 20/E0/1000.0
  +Ef = 10/E0/1000.0
   QSYS_DEBUG_LEVEL = 2 ! show more info
   
   
  
positive.code: !    enddo
  !    enddo
  
      call gemm(invGamma,tmpSkm,dSkm)
      print"(A,6f12.4)","G(p,p)=",G0(p,p),G0(p,p)/ldos
  
  
      waveN = 0
      do m = 1 , qt%leads(1)%no_in_modes
          do i = 1 , qt%leads(2)%no_sites
              la = qt%leads(2)%l2g(i,1)
              ls = qt%leads(2)%l2g(i,2)
              lg = qt%qsystem%atoms(la)%globalIDs(ls);
              waveN(m,lg) = dSkm(i,m)
  !            waveN(m,lg) = tmpSkm(i,m)
          enddo
      enddo
      waveN = waveN + waveL
  
      total_Tn = 0
      total_Rn = 0
      qt%Tn    = 0
      qt%Rn    = 0
      do m = 1 , qt%leads(1)%no_in_modes
  
      do i = 1 , qt%no_leads
          qt%leads(i)%totalT = 0
          ! Skip tranmission from input lead
          if( i /= leadID) then
              call qt%leads(i)%calculate_Tnm(qt%qsystem%atoms,waveN(m,:))
              total_Tn = total_Tn  + qt%leads(i)%modeT
          else
              call qt%leads(i)%calculate_Tnm(qt%qsystem%atoms,waveN(m,:),m)
              total_Rn = total_Rn  + qt%leads(i)%modeT
          endif
          qt%leads(i)%totalT = qt%leads(i)%totalT + qt%leads(i)%modeT
      enddo
          qt%Tn(m) = total_Tn
          qt%Rn(m) = total_Rn
      enddo
      write(112,"(100e20.6)"),ytip*dx,total_Tn
      print"(100e20.6)",ytip*dx,total_Tn
  
      ! calculate QdagUdag * dS
      no_modes = qt%leads(2)%no_in_modes
  
negative.code: 
  QSYS_DEBUG_LEVEL = 0
  open(unit=111,file=output_folder//"sgmT.dat")
  open(unit=112,file=output_folder//"sgmP1.dat")
  xtip = nx/2-5
  do ytip = 1 , ny
  do xtip = nx/2-nx/4 , nx/2+nx/4
      Vpot = 0
      Vtip = Ef/1.0
  !    T1 = get_clock()
      Vpot(gindex(xtip,ytip)) = Vtip
  !    call qt%qsystem%update_lattice(c_matrix=coupling)
  !    call qt%calculate_modes(Ef)
  !    call qt%solve(1,Ef)
  !    write(111,"(100e20.6)"),ytip*dx,sum(qt%Tn(:))
  !    print*,"time full:",get_clock()-T1
      T1 = get_clock()
      call calc_sgm()
  !    print*,"time fast:",get_clock()-T1
  enddo
      write(112,*),""
  enddo
  close(111)
  
  