From 59940c18263a854a92dac57a958c42ab6702d793 Mon Sep 17 00:00:00 2001
From: smogunov <smogunov@c92efa57-630b-4861-b058-cf58834340f0>
Date: Wed, 9 Mar 2005 17:44:51 +0000
Subject: [PATCH] Check on the transmission matrix unitarity is added
 (A.Smogunov)

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1697 c92efa57-630b-4861-b058-cf58834340f0
---
 PWCOND/Makefile     |  1 +
 PWCOND/sunitary.f90 | 48 +++++++++++++++++++++++++++++++++++++++++++++
 PWCOND/transmit.f90 | 23 +++++++++++++++++++++-
 3 files changed, 71 insertions(+), 1 deletion(-)
 create mode 100644 PWCOND/sunitary.f90

diff --git a/PWCOND/Makefile b/PWCOND/Makefile
index 856b92aa1..b53a75f36 100644
--- a/PWCOND/Makefile
+++ b/PWCOND/Makefile
@@ -42,6 +42,7 @@ scatter_back.o \
 scatter_forw.o \
 summary_band.o \
 summary_tran.o \
+sunitary.o \
 transmit.o
 
 MODULES = \
diff --git a/PWCOND/sunitary.f90 b/PWCOND/sunitary.f90
new file mode 100644
index 000000000..0b98d1127
--- /dev/null
+++ b/PWCOND/sunitary.f90
@@ -0,0 +1,48 @@
+!
+! Copyright (C) 2003 A. Smogunov
+! This file is distributed under the terms of the
+! GNU General Public License. See the file `License'
+! in the root directory of the present distribution,
+! or http://www.gnu.org/copyleft/gpl.txt .
+!
+subroutine sunitary(nchanl, nchanr, smat, info)
+!
+! It performs a check of unitarity of the scattering matrix
+! S = {R_ij, T_ij}
+! \sum R_ki^* R_kj + \sum T_ki^* R_kj = \delta_ij
+!
+  USE kinds, ONLY : DP
+  USE io_global,  ONLY :  stdout
+implicit none
+
+  integer ::  info, i, j, k
+  integer ::  nchanl, &  ! number of channels in the left tip
+              nchanr     !     ------------         right tip
+  real(kind=DP), parameter :: eps=1.d-4
+  real(kind=DP) :: raux, raux1
+  complex(kind=DP) ::   &
+              s, smat(nchanl+nchanr, nchanl)    ! S matrix
+
+  info = 0
+  do i = 1, nchanl
+    do j = 1, nchanl
+      s = 0.d0
+      do k = 1, nchanl
+       s = s + conjg(smat(k,i))*smat(k,j)
+      enddo
+      do k = 1, nchanr
+       s = s + conjg(smat(nchanl+k,i))*smat(nchanl+k,j)
+      enddo
+      raux = sqrt(real(s)**2 + aimag(s)**2)
+      raux1 = raux
+      if(i.eq.j) raux1 = abs(raux1-1.d0)
+      if(raux1.gt.eps) then
+        write(stdout, '(2i5,f12.7)') i, j, raux
+        info = 1
+      endif
+    enddo
+  enddo
+
+  return
+end subroutine sunitary
+
diff --git a/PWCOND/transmit.f90 b/PWCOND/transmit.f90
index ef6ee269b..b3783266a 100644
--- a/PWCOND/transmit.f90
+++ b/PWCOND/transmit.f90
@@ -30,7 +30,7 @@ implicit none
   complex(kind=DP) :: x1, x2
   complex(kind=DP), allocatable :: amat(:,:), vec1(:,:), &
                      tmat(:,:), veceig(:,:), zps_nc(:,:), &
-                     vec2(:,:)
+                     vec2(:,:), smat(:,:)
 
   eev = earr(ien)
   ntran=4*n2d+npol*(norbs+nocrosl+nocrosr)
@@ -221,6 +221,27 @@ implicit none
     enddo
     WRITE( stdout,'(15x,f9.5)') tj
   enddo
+
+!
+! Check for S matrix unitarity
+!
+!
+! elements of S_ij = (r_ij, t_ij) and T_ij = (t_ij) matrices
+!
+  allocate( smat(nchanl+nchanr, nchanl) )
+  do n=1, nchanl
+    do ig=1, nchanl
+      smat(ig,n) = vec1(2*n2d+npol*norbs+ig,n)
+    enddo
+    do ig=1, nchanr
+      smat(nchanl+ig,n) = tmat(ig,n)
+    enddo
+  enddo
+  call sunitary(nchanl, nchanr, smat, info)
+  call errore('transmit','S matrix is not unitary', &
+                                     -abs(info))
+  deallocate( smat )
+
 !
 ! eigenchannel decomposition
 !