Experimental: pwcond calculates the scattering and Bloch states on a 3D-FFT

mesh and writes them in a file readable by the pp programs. git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3490 c92efa57-630b-4861-b058-cf58834340f0
2006-10-25 15:35:53 +00:00 · 2006-10-25 15:35:53 +00:00 · b6afc58105
parent 0b1d3aeb5a
commit b6afc58105
6 changed files with 152 additions and 79 deletions
--- a/PWCOND/Makefile
+++ b/PWCOND/Makefile
@ -4,7 +4,7 @@ include ../make.sys

 PPOBJS = \
 ../PP/start_postproc.o \
-../PP/xsf.o \
+../PP/plot_io.o \
 ../PP/stop_pp.o

 PWCONDOBJS = \
@ -39,7 +39,8 @@ scatter_forw.o \
 summary_band.o \
 summary_tran.o \
 sunitary.o \
-transmit.o
+transmit.o \
+write_states.o

 MODULES = \
 ../Modules/atom.o \
--- a/PWCOND/allocate_cond.f90
+++ b/PWCOND/allocate_cond.f90
@ -33,10 +33,7 @@ subroutine allocate_cond
  allocate( fund0(n2d, 2*n2d) )
  allocate( fund1(n2d, 2*n2d) ) 

-  IF (lorb) THEN
-     allocate( funz0(n2d, 2*n2d+norbf*npol, nrzpl) )
-     allocate( kfunz(n2d, 2*(n2d+npol*nocrosl), nrzpl) )
-  END IF
+  IF (lorb) allocate( funz0(n2d, 2*n2d+norbf*npol, MAX(nrzpl,nrzps)) )

  IF (norbf>0) THEN
     allocate( funl0(n2d, norbf*npol) )     
--- a/PWCOND/compbs.f90
+++ b/PWCOND/compbs.f90
@ -32,7 +32,7 @@ subroutine compbs(lleft, nocros, norb, nchan, kval, kfun,  &
     lleft       ! 1/0 if it is left/right tip    
  integer :: ik, i, j, kin, kfin, ig, info, lam, n, iorb, iorb1, &
             iorb2, aorb, borb, nt, startk, lastk, nchan, nb, ih, &
-             ih1, ij, is, js, ichan 
+             ih1, ij, is, js, ichan
  REAL(DP), PARAMETER :: eps=1.d-8
  REAL(DP) :: raux, DDOT
  REAL(DP), ALLOCATABLE :: zpseu(:,:,:), zps(:,:)
@ -354,73 +354,14 @@ subroutine compbs(lleft, nocros, norb, nchan, kval, kfun,  &

  endif

-IF (lorb) THEN
-!
-!    Calculate the Bloch functions in all z points
-!
-   raux=0.d0
-   DO ik=1,nrzpl
-      CALL ZGEMM('n', 'n', n2d, nchan, 2*n2d+npol*norb, one, funz0(1,1,ik), &
-                   n2d, vec, 2*n2d+npol*norb, zero, kfunz(1,1,ik), n2d)
-      raux=raux+DDOT(2*n2d,kfunz(1,1,ik),1,kfunz(1,1,ik),1)
-   END DO
-   raux=raux*omega/nrzpl
-   raux=1.d0/sqrt(raux)
-   DO ik=1,nrzpl
-      call DSCAL(2*n2d,raux,kfunz(1,1,ik),1)
-   ENDDO
-!
-!    Transform and save all propagating functions. One function after the other.
-!
-   DO ik=1,nrzpl   
-      DO ichan=1,nchan
-         fung=(0.d0,0.d0)
-         DO ig=1,npol*ngper
-            DO mu=1,n2d
-               fung(ig)=fung(ig) + kfunz(mu,ichan,ik) * newbg(ig,mu)
-            END DO
-         END DO
-
-         funr=(0.d0,0.d0)
-         DO ig=1,ngper*npol
-            funr(nl_2d(ig))=fung(ig)
-         END DO
-
-         CALL cft3(funr,nrx,nry,1,nrx,nry,1,1)
-!
-!  Copy the wavefunction for printing, using the boundary conditions
-!
-         DO ix=1,nrx+1
-            iix=ix
-            IF (ix>nrx) iix=1
-            DO jx=1,nry+1
-               jjx=jx
-               IF (jx>nry) jjx=1
-               ig=ix+(jx-1)*(nrx+1)
-               iig=iix+(jjx-1)*nrx
-               auxfr(ig)=REAL(funr(iig))**2+AIMAG(funr(iig))**2
-            END DO
-         END DO
-!  
-         ounit = 34
-         IF (ik>9) THEN
-            write(filename,'("wfc_z.",i2)') ik
-         ELSE
-            write(filename,'("wfc_z.",i1)') ik
-         ENDIF
-         OPEN (UNIT=ounit, FILE=filename, FORM='formatted', STATUS='unknown')
-
-         x0=0.d0
-         CALL xsf_struct (alat, at, 1, x0, 'Al1', 1, ounit)
-         CALL xsf_datagrid_2d (auxfr,nrx+1,nry+1,1.d0,1.d0,x0,at(1,1),at(1,2),&
-                                                alat,ounit)
-!         do ix=1,nrx+1
-!            write(ounit,*) ix, REAL(auxfr(ix))
-!         ENDDO
-         CLOSE(ounit)
-      END DO
-   END DO
-END IF
+  IF (lorb.and.ikind==0) THEN
+  !
+  !    Calculate the Bloch functions in all z points
+  !
+     ounit=34
+     CALL write_states(nrzpl, n2d, norb, norbf, nchan, &
+                                   nrx, nry, ounit, funz0, vec, .TRUE.)
+  END IF

  deallocate(amat)
  deallocate(bmat)
--- a/PWCOND/condcom.f90
+++ b/PWCOND/condcom.f90
@ -196,7 +196,6 @@ MODULE cb_cond
  COMPLEX(DP), ALLOCATABLE :: &
       kvall(:),       &!  k           for the left lead 
       kfunl(:,:),     &!  phi_k(z=d)  for the left lead
-       kfunz(:,:,:),   &!  phi_k(z)  for all slabs
       kfundl(:,:),    &!  phi_k'(z=d) for the left lead
       kintl(:,:),     &!  integral of phi_k with beta-fun.
       kcoefl(:,:),    &!  coeff. of phi_k over nonloc. fun.
--- a/PWCOND/transmit.f90
+++ b/PWCOND/transmit.f90
@ -23,14 +23,14 @@ subroutine transmit(ik, ien)
 implicit none

  integer :: ik, ien, n, iorb, iorb1, iorb2, iorba, ipol, nt, &
-             ih, ih1, ig, ntran, ij, is, js, info
+             ih, ih1, ig, ntran, ij, is, js, ichan, ounit, info
  integer, allocatable :: ipiv(:)
  real(DP) :: tk, tj, tij, eev
  real(DP), allocatable :: zps(:,:), eigen(:)
  complex(DP) :: x1, x2, xi1(2)
  complex(DP), allocatable :: amat(:,:), vec1(:,:), &
                     tmat(:,:), veceig(:,:), zps_nc(:,:), &
-                     vec2(:,:), smat(:,:)
+                     vec2(:,:), smat(:,:), vec(:,:)

  eev = earr(ien)
  ntran=4*n2d+npol*(norbs+nocrosl+nocrosr)
@ -303,8 +303,9 @@ implicit none
               xi1(ipol) = xi1(ipol)+intw2(iorba, ig)*vec2(2*n2d+ig, n)
            enddo
         enddo
-         write(stdout,'(2i5,2f20.12)') n, iorb-orbj_in+1,   &
-                        ( DBLE(xi1(ipol))**2+AIMAG(xi1(ipol))**2,ipol=1,npol)
+         write(stdout,'(2i5,2f16.12,2x,2f16.12)') n, iorb-orbj_in+1,   &
+                                      (xi1(ipol),ipol=1,npol)
+
      enddo
     endif
    enddo
@ -313,6 +314,18 @@ implicit none

  endif

+  IF (lorb) THEN
+     ALLOCATE(vec(2*n2d+npol*norbs,nchanl))
+     DO ichan=1,nchanl
+        DO ig=1, 2*n2d+npol*norbs
+           vec(ig,ichan)=vec1(ig,ichan)
+        END DO
+     END DO
+     ounit=34
+     CALL write_states(nrzps,n2d,norbs,norbf,nchanl,nrx,nry, &
+                                       ounit,funz0,vec,.FALSE.)
+     DEALLOCATE(vec)
+  END IF

  deallocate(ipiv)
  deallocate(amat)
--- a/PWCOND/write_states.f90
+++ b/PWCOND/write_states.f90
@ -0,0 +1,122 @@
+!
+! Copyright (C) 2006 QUANTUM-ESPRESSO
+! 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 write_states(nrzp,n2d,norb,norbf,nchan,nrx,nry,ounit, &
+                                                  funz0,vec,norm_flag)
+!
+! Using the basis functions obtained by scatter_forw (funz0), and
+! the coefficients obtained in compbs or in transmit (vec) it computes
+! the states on a three dimensional grid and writes them in files
+! readable by the post-processing programs. 
+! If norm_flag is .true. the states are normalized in such a way
+! that their square modulus integrated over all volume is one.
+!
+#include "f_defs.h"
+ USE kinds, ONLY : DP
+ USE noncollin_module, ONLY : npol
+ USE cond,      ONLY : ngper, newbg, nl_2d
+ USE pwcom,     ONLY : omega, ibrav, celldm, gcutm, dual, ecutwfc, at
+ USE ions_base, ONLY : ityp, zv, tau, atm
+
+ IMPLICIT NONE
+ INTEGER :: nrzp,n2d,norb,norbf,nchan,nrx,nry,ounit 
+ COMPLEX(DP) :: funz0(n2d, 2*n2d+norbf*npol, nrzp), vec(2*n2d+npol*norb,nchan)
+ LOGICAL :: norm_flag
+ INTEGER :: ik, ig, ir, mu, ichan, ios
+ REAL(DP) :: DDOT
+ COMPLEX(DP), PARAMETER :: one=(1.d0,0.d0), zero=(0.d0,0.d0)
+ CHARACTER(LEN=50) :: filename
+ CHARACTER(LEN=75) :: title_here
+
+ INTEGER :: ix, jx
+ COMPLEX(DP), ALLOCATABLE :: funr(:), fung(:), kfunz(:,:,:)
+
+ REAL(DP), ALLOCATABLE :: aux_plot(:), norm(:)
+!
+!    Calculate the functions in all z points
+!
+ ALLOCATE(funr(nrx*nry))
+ ALLOCATE(fung(npol*ngper))
+ ALLOCATE(kfunz(n2d,nchan,nrzp))
+ DO ik=1,nrzp
+    CALL ZGEMM('n', 'n', n2d, nchan, 2*n2d+npol*norb, one, funz0(1,1,ik), &
+                n2d, vec, 2*n2d+npol*norb, zero, kfunz(1,1,ik), n2d)
+ END DO
+ IF (norm_flag) THEN
+    ALLOCATE (norm(nchan))
+    norm=0.d0
+    DO ichan=1,nchan
+       DO ik=1,nrzp
+          norm(ichan) = norm(ichan) + &
+                     DDOT(2*n2d,kfunz(1,ichan,ik),1,kfunz(1,ichan,ik),1)
+       END DO
+    END DO
+    norm=norm*omega/nrzp
+    norm=1.d0/SQRT(norm)
+    DO ik=1,nrzp
+       DO ichan=1,nchan
+          CALL DSCAL(2*n2d,norm(ichan),kfunz(1,ichan,ik),1)
+       END DO 
+    END DO
+    DEALLOCATE (norm)
+ END IF
+!
+!    Transform and save all propagating functions. One function after the other.
+!
+ ALLOCATE(aux_plot(nrx*nry*nrzp))
+ DO ichan=1,nchan
+    ounit = 34
+    IF (ichan>9) THEN
+      write(filename,'("wfc_n.",i2)') ichan
+    ELSE
+      write(filename,'("wfc_n.",i1)') ichan
+    ENDIF
+    filename=TRIM(filename)
+    OPEN (UNIT=ounit, FILE=filename, FORM='formatted', &
+                           STATUS='unknown', ERR=100, IOSTAT=ios)
+100 CALL errore('write_states','opening file'//filename,ABS(ios))
+
+    aux_plot=(0.d0,0.d0)
+    DO ik=1,nrzp
+       fung=(0.d0,0.d0)
+       DO ig=1,npol*ngper
+          DO mu=1,n2d
+             fung(ig)=fung(ig) + kfunz(mu,ichan,ik) * newbg(ig,mu)
+          END DO
+       END DO
+
+       funr=(0.d0,0.d0)
+       DO ig=1,ngper*npol
+          funr(nl_2d(ig))=fung(ig)
+       END DO
+
+       CALL cft3(funr,nrx,nry,1,nrx,nry,1,1)
+
+       DO ix=1,nrx
+          DO jx=1,nry
+             ir=ix + (jx - 1) * nrx + (ik - 1) * nrx * nry
+             ig=ix+(jx-1)*nrx
+             aux_plot(ir)=REAL(funr(ig))**2+AIMAG(funr(ig))**2
+          END DO
+       END DO
+    END DO
+    CLOSE(ounit)
+
+    title_here='written by pwcond'
+
+    CALL plot_io (TRIM(filename), title_here, nrx, nry, nrzp,       &
+                  nrx, nry, nrzp, 0, 0, ibrav, celldm, at, gcutm,   &
+                  dual, ecutwfc, 7, atm, ityp, zv, tau, aux_plot, + 1)
+ END DO
+ DEALLOCATE(aux_plot)
+ DEALLOCATE(kfunz)
+ DEALLOCATE(funr)
+ DEALLOCATE(fung)
+
+ RETURN
+END SUBROUTINE write_states