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: 
  
  tmp2GU = -transpose(qt%leads(2)%modes(M_OUT,:,:))
  call gemm(qt%leads(2)%GammaVecs,tmp2GU,tmpGU)
  
  
  
  call zalg_SVD_QL(tmpGU,L,QdagUdag)
  !print*,"cond(L)       =",alg_cond(L)
  print*,"cond(L(modes))=",alg_cond(L(1:qt%leads(2)%no_in_modes,1:qt%leads(2)%no_in_modes))
  !print*,"cond(QdagUdag)=",alg_cond(QdagUdag)
  !
  !call gemm(QdagUdag,L,tmp2GU,transa="C")
  !tmpGU = transpose(qt%leads(2)%modes(M_OUT,:,:))
  !call inverse_matrix(qt%leads(2)%no_sites,tmpGU)
  !call gemm(tmp2GU,tmpGU,QdagUdag)
  !print*,"sum(GammaVec2)=",sum(QdagUdag)
  !stop
  !call write_matrix(L,output_folder//"L.dat","(500f10.4)")
  
  
  
  
  print*,"T=",sum(abs(smatrix0(1,2)%Tnm)**2),"R=",sum(abs(smatrix0(1,1)%Tnm)**2)
  
negative.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
  