From b6afc58105c38626a468984e0260fc4f040de8a3 Mon Sep 17 00:00:00 2001 From: dalcorso Date: Wed, 25 Oct 2006 15:35:53 +0000 Subject: [PATCH] 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 --- PWCOND/Makefile | 5 +- PWCOND/allocate_cond.f90 | 5 +- PWCOND/compbs.f90 | 77 +++--------------------- PWCOND/condcom.f90 | 1 - PWCOND/transmit.f90 | 21 +++++-- PWCOND/write_states.f90 | 122 +++++++++++++++++++++++++++++++++++++++ 6 files changed, 152 insertions(+), 79 deletions(-) create mode 100644 PWCOND/write_states.f90 diff --git a/PWCOND/Makefile b/PWCOND/Makefile index 1e36bb1ed..35186ac3d 100644 --- 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 \ diff --git a/PWCOND/allocate_cond.f90 b/PWCOND/allocate_cond.f90 index 24f915fd0..3cfec4f99 100644 --- 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) ) diff --git a/PWCOND/compbs.f90 b/PWCOND/compbs.f90 index a1aa8b6ef..55dfb7f9d 100644 --- 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) diff --git a/PWCOND/condcom.f90 b/PWCOND/condcom.f90 index 7152675e7..56947d224 100644 --- 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. diff --git a/PWCOND/transmit.f90 b/PWCOND/transmit.f90 index 1c8df61a5..d57a10fc7 100644 --- 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) diff --git a/PWCOND/write_states.f90 b/PWCOND/write_states.f90 new file mode 100644 index 000000000..c48088512 --- /dev/null +++ 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