scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Modified the transmission code PWCOND so that now it works 

with the separate files describing the leads and the scattering region (more accurate
description of the complex bands of the leads) as well as with the unique file
containing the scattering region + the leads. Added an option allowing to save on
a text file (machine independent) the data needed for PWCOND run.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1590 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
smogunov 2005-01-28 17:10:20 +00:00
parent f50f5fdda2
commit 064d48a81d
28 changed files with 2156 additions and 1594 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
PWCOND/Makefile
View File

@ -9,13 +9,13 @@ PPOBJS = \
PWCONDOBJS = \
allocate_cond.o \
allocate_cond_2.o \
bessj.o \
cft3sp.o \
compbs.o \
compbs_2.o \
condcom.o \
condmain.o \
cond_out.o \
do_cond.o \
eigenchnl.o \
form_zk.o \
@ -27,6 +27,7 @@ gep_x.o \
gramsh.o \
hev_ab.o \
init_cond.o \
init_orbitals.o \
init_gper.o \
integrals.o \
jbloch.o \
@ -37,18 +38,17 @@ poten.o \
print_clock_pwcond.o \
rotate.o \
rotproc.o \
save_cond.o \
scatter_back.o \
scatter_forw.o \
slabcpu.o \
summary_band.o \
summary_tran.o \
transmit.o
MODULES = \
../Modules/atom.o \
../Modules/basic_algebra_routines.o \
../Modules/berry_phase.o \
../Modules/cell_base.o \
../Modules/clocks.o \
../Modules/check_stop.o \
@ -119,8 +119,10 @@ PWOBJS = \
../PW/checksym.o \
../PW/cinitcgg.o \
../PW/cinitcgg_nc.o \
../PW/clean_pw.o \
../PW/compute_dip.o \
../PW/compute_rho.o \
../PW/constraints_module.o \
../PW/coset.o \
../PW/cryst_to_car.o \
../PW/cubicsym.o \
@ -193,6 +195,7 @@ PWOBJS = \
../PW/openfil.o \
../PW/ortho.o \
../PW/orthoatwfc.o \
../PW/output_tau.o \
../PW/para.o \
../PW/potinit.o \
../PW/print_clock_pw.o \

250
PWCOND/allocate_cond.f90
View File

@ -1,229 +1,59 @@
!
! Copyright (C) 2003 A. Smogunov
! 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 .
!
! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC).
!
subroutine allocate_cond
!
! This subroutine allocates the variables for PWCOND
! This subroutine allocates some needed variables
!
#include "f_defs.h"
USE ions_base, ONLY : nat, ityp, ntyp => nsp, tau
use pwcom
USE noncollin_module, ONLY : noncolin, npol
USE uspp_param, ONLY : nbrx, nbeta, lll, betar, tvanp
use atom, only: mesh, r
USE noncollin_module, ONLY : npol
use cond
implicit none
integer :: mmax, nt, k, iw, ib, ir, na, naux
real(kind=DP), parameter :: eps=1.d-8
real(kind=DP) :: &
bmax, & ! maximum of beta function
ledge, & ! left edge of the orbital
redge, & ! right edge of the orbital
dz1, dz2, & ! auxiliary
epsbeta
real(kind=DP), allocatable :: dwid(:)
!
! FFT parameters for local potential
!
nrx=nr1s
nry=nr2s
nrz=nr3s+dnslab
if(nrz/2*2.eq.nrz) nrz=nrz+1
allocate( newbg(ngper*npol, n2d) )
allocate( psiperl( n2d, n2d, nrzpl ) )
allocate( zkl( n2d, nrzpl ) )
allocate( zkrl( n2d, nrzpl ) )
allocate( psipers( n2d, n2d, nrzps ) )
allocate( zks( n2d, nrzps ) )
allocate( zkrs( n2d, nrzps ) )
allocate( psiperr( n2d, n2d, nrzpr ) )
allocate( zkr( n2d, nrzpr ) )
allocate( zkrr( n2d, nrzpr ) )
allocate ( z(nrz+1) )
allocate ( rsph(nbrx, ntypx) )
allocate ( dwid(6) )
!
! slab parameters
!
zl=at(3,3)
sarea=abs(at(1,1)*at(2,2)-at(2,1)*at(1,2))*alat**2
if (abs(ecut2d).le.eps) ecut2d=ecutwfc
allocate( fun0(n2d, 2*n2d) )
allocate( fun1(n2d, 2*n2d) )
allocate( fund0(n2d, 2*n2d) )
allocate( fund1(n2d, 2*n2d) )
if (ikind.gt.2.or.ikind.le.-1) &
call errore ('init_cond','ikind is not in allowed range',1)
if (ikind.eq.0) then
bdr1=bdl2
bdr2=bdl2
endif
if (ikind.eq.1) then
bdr1 = zl
bdr2 = zl
allocate( funl0(n2d, norbf*npol) )
allocate( funl1(n2d, norbf*npol) )
allocate( fundl0(n2d, norbf*npol) )
allocate( fundl1(n2d, norbf*npol) )
allocate( intw1(norbf*npol, 2*n2d) )
allocate( intw2(norbf*npol, norbf*npol) )
allocate( kvall(2*(n2d+npol*nocrosl)) )
allocate( kfunl(n2d, 2*(n2d+npol*nocrosl)) )
allocate( kfundl(n2d, 2*(n2d+npol*nocrosl)) )
allocate( kintl(nocrosl*npol, 2*(n2d+npol*nocrosl)) )
allocate( kcoefl(nocrosl*npol, 2*(n2d+npol*nocrosl)) )
if(ikind.ne.0) then
allocate( kvalr(2*(n2d+npol*nocrosr)) )
allocate( kfunr(n2d, 2*(n2d+npol*nocrosr)) )
allocate( kfundr(n2d, 2*(n2d+npol*nocrosr)) )
allocate( kintr(nocrosr*npol, 2*(n2d+npol*nocrosr)) )
allocate( kcoefr(nocrosr*npol, 2*(n2d+npol*nocrosr)) )
endif
dz1=zl/nrz
dwid(1)=0.d0
dwid(2)=bdl1
dwid(3)=bdl2
dwid(4)=bdr1
dwid(5)=bdr2
dwid(6)=zl
nrz=0
do iw=2, 6
naux=(dwid(iw)-dwid(1)+dz1*0.5d0)/dz1-nrz
if (naux.gt.0) then
dz2=(dwid(iw)-dwid(iw-1))/naux
do k=1, naux
z(nrz+k)=dwid(iw-1)+(k-1)*dz2
enddo
nrz=nrz+naux
endif
enddo
z(nrz+1)=zl
!
! to determine radii of nonlocal spheres
!
epsbeta=1.d-4
mmax = 0
do nt=1, ntyp
do ib=1, nbeta(nt)
mmax = max(mmax, lll(ib, nt))
bmax=0.d0
do ir=2, mesh(nt)
bmax=max(bmax, betar(ir,ib,nt)/r(ir,nt))
enddo
ir=mesh(nt)
do while (abs(betar(ir,ib,nt)/r(ir,nt)).le.epsbeta*bmax)
ir=ir-1
enddo
rsph(ib,nt)=r(ir,nt)/alat
enddo
enddo
if (mmax.gt.3) call errore ('allocate','for l>3 -orbitals &
& the method is not implemented',1)
!
! We set up the radii of US spheres to be the same (to avoid
! the problem with the spheres crossing or not the boundaries)
!
do nt=1, ntyp
if (tvanp(nt)) then
bmax=0.d0
do ib=1, nbeta(nt)
bmax=max(bmax, rsph(ib,nt))
enddo
do ib=1, nbeta(nt)
rsph(ib,nt)=bmax
enddo
endif
enddo
!
! Move all atoms into the unit cell
!
do na=1, nat
tau(3,na) = tau(3,na)/zl
tau(3,na) = tau(3,na) - int(tau(3,na))
if(tau(3,na).le.eps) tau(3,na)=tau(3,na)+1.d0
tau(3,na) = tau(3,na)*zl
enddo
!
! Compute the number of orbitals inside the cell and
! crossing the boundary
!
noinsl=0
nocrosl=0
noinss=0
noinsr=0
nocrosr=0
norbs=0
!
! for left tip
!
do na=1, nat
nt=ityp(na)
do ib=1, nbeta(nt)
ledge=tau(3,na)-rsph(ib,nt)
redge=tau(3,na)+rsph(ib,nt)
if (ledge.le.bdl1.and.redge.gt.bdl2) &
call errore ('init_cond','The unit cell in left',1)
if (tau(3,na).gt.bdl1+eps.and.tau(3,na).le.bdl2+eps) then
if (ledge.gt.bdl1.and.redge.le.bdl2) then
noinsl=noinsl+2*lll(ib,nt)+1
else
nocrosl=nocrosl+2*lll(ib,nt)+1
endif
endif
enddo
enddo
norb=2*nocrosl+noinsl
norbf=norb
!
! If scattering problem
!
if (ikind.ne.0) then
!
! Scatt. region
!
do na=1, nat
nt=ityp(na)
do ib=1, nbeta(nt)
ledge=tau(3,na)-rsph(ib,nt)
redge=tau(3,na)+rsph(ib,nt)
if (ledge.le.bdl2.and.redge.gt.bdr1) &
call errore ('init_cond','The unit cell in scatt.',1)
if (ledge.gt.bdl2.and.redge.le.bdr1) &
noinss=noinss+2*lll(ib,nt)+1
enddo
enddo
norb=norb+noinss
!
! Right tip
!
if (ikind.eq.1) then
!
! If the tips are equal
!
nocrosr=nocrosl
noinsr=noinsl
norb=norb+nocrosr
!
! If the tips are different
!
else
do na=1, nat
nt=ityp(na)
do ib=1, nbeta(nt)
ledge=tau(3,na)-rsph(ib,nt)
redge=tau(3,na)+rsph(ib,nt)
if (ledge.le.bdr1.and.redge.gt.bdr2) &
call errore ('init_cond','The unit cell in right',1)
if (ledge.gt.bdr1.and.redge.le.bdr2) &
noinsr=noinsr+2*lll(ib,nt)+1
if (ledge.le.bdr1.and.redge.gt.bdr1) &
nocrosr=nocrosr+2*lll(ib,nt)+1
enddo
enddo
norb=norb+2*nocrosr+noinsr
norbf=max(norbf, 2*nocrosr+noinsr)
endif
norbs=nocrosl+noinss+nocrosr
norbf=max(norbf, norbs)
endif
allocate( vppot(nrz, nrx * nry, npol, npol) )
allocate( itnew(norb) )
allocate( nbnew(norb) )
allocate( natih(norb, 2) )
allocate( ls(norb) )
allocate( mnew(norb) )
allocate( cross(norb, nrz) )
allocate( taunew(3, norb) )
allocate( zpseu(norb, norb, nspin) )
if (noncolin) allocate(zpseu_nc(norb, norb, nspin))
allocate( nkofz(nrz) )
deallocate(dwid)
return
return
end subroutine allocate_cond

53
PWCOND/allocate_cond_2.f90
View File

@ -1,53 +0,0 @@
! 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 .
!
! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC).
!
subroutine allocate_cond_2
!
! This subroutine allocates the remaining variables
! after reduction (ngper --> n2d) of XY basis set
!
#include "f_defs.h"
USE noncollin_module, ONLY : npol
use cond
implicit none
allocate( newbg(ngper*npol, n2d) )
allocate( psiper( n2d, n2d, nrzp ) )
allocate( zk( n2d, nrzp ) )
allocate( zkr( n2d, nrzp ) )
allocate( fun0(n2d, 2*n2d) )
allocate( fun1(n2d, 2*n2d) )
allocate( fund0(n2d, 2*n2d) )
allocate( fund1(n2d, 2*n2d) )
allocate( funl0(n2d, norbf*npol) )
allocate( funl1(n2d, norbf*npol) )
allocate( fundl0(n2d, norbf*npol) )
allocate( fundl1(n2d, norbf*npol) )
allocate( intw1(norbf*npol, 2*n2d) )
allocate( intw2(norbf*npol, norbf*npol) )
allocate( kvall(2*(n2d+npol*nocrosl)) )
allocate( kfunl(n2d, 2*(n2d+npol*nocrosl)) )
allocate( kfundl(n2d, 2*(n2d+npol*nocrosl)) )
allocate( kintl(nocrosl*npol, 2*(n2d+npol*nocrosl)) )
allocate( kcoefl(nocrosl*npol, 2*(n2d+npol*nocrosl)) )
if(ikind.ne.0) then
allocate( kvalr(2*(n2d+npol*nocrosr)) )
allocate( kfunr(n2d, 2*(n2d+npol*nocrosr)) )
allocate( kfundr(n2d, 2*(n2d+npol*nocrosr)) )
allocate( kintr(nocrosr*npol, 2*(n2d+npol*nocrosr)) )
allocate( kcoefr(nocrosr*npol, 2*(n2d+npol*nocrosr)) )
endif
return
end subroutine allocate_cond_2

236
PWCOND/compbs.f90
View File

@ -1,5 +1,3 @@
!
! Copyright (C) 2003 A. Smogunov
! This file is distributed under the terms of the
@ -10,32 +8,32 @@
! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC).
!
!
subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
nchan, kval, kfun, kfund, kint, kcoef)
subroutine compbs(lleft, nocros, norb, nchan, kval, kfun, &
kfund, kint, kcoef)
!
! Using the basis functions obtained by scatter_forw it computes
! the complex band structure (CBS) of the tip ( zin<z<zfin ).
! the complex band structure (CBS) of the lead.
! Some variables needed for wave-function matching in transmission
! calculation are constructed and saved.
!
#include "f_defs.h"
use pwcom
USE noncollin_module, ONLY : noncolin, npol
USE uspp_param, ONLY: tvanp
use spin_orb, only : lspinorb
use lsda_mod, only: nspin
use cond
implicit none
integer :: &
nocros, & ! number of orbitals crossing the boundary
noins, & ! number of interior orbitals
norbnow, & ! total number of orbitals
orbin, & ! number of 1-st orbital in full list
lright ! 1/0 if it is right/left tip
norb, & ! total number of orbitals
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
real(kind=DP), parameter :: eps=1.d-8
real(kind=DP) :: raux, zin, zfin, dz, DDOT
real(kind=DP), allocatable :: zps(:,:)
real(kind=DP) :: raux, DDOT
real(kind=DP), allocatable :: zpseu(:,:,:), zps(:,:)
complex(kind=DP), parameter :: cim=(0.d0,1.d0)
complex(kind=DP) :: x1, x2, x3, x4, y1, y2, y3, y4, &
kval(2*(n2d+npol*nocros)), kfun(n2d,2*(n2d+npol*nocros)), &
@ -43,92 +41,75 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
kcoef(nocros*npol,2*(n2d+npol*nocros)), &
kfund(n2d,2*(n2d+npol*nocros))
complex(kind=DP), allocatable :: amat(:,:), bmat(:,:), vec(:,:), &
zps_nc(:,:), aux(:,:)
zpseu_nc(:,:,:,:), &
zps_nc(:,:), aux(:,:)
complex(kind=DP), parameter :: one=(1.d0,0.d0), zero=(0.d0,0.d0)
LOGICAL :: exst
call start_clock('compbs')
noins = norbnow-2*nocros
allocate( amat( (2*n2d+npol*norbnow), (2*n2d+npol*norbnow) ) )
allocate( bmat( (2*n2d+npol*norbnow), (2*n2d+npol*norbnow) ) )
allocate( vec( (2*n2d+npol*norbnow), 2*(n2d+npol*nocros) ) )
allocate( zps( norbnow, norbnow ) )
allocate( aux( n2d, 2*n2d+npol*norbnow))
if (noncolin) allocate( zps_nc( norbnow*npol, norbnow*npol ) )
noins = norb-2*nocros
if(lleft.eq.1) then
if (noncolin) then
allocate( zpseu_nc(2,norb,norb,nspin) )
zpseu_nc = zpseul_nc
else
allocate( zpseu(2,norb,norb) )
zpseu = zpseul
endif
else
if (noncolin) then
allocate( zpseu_nc(2,norb,norb,nspin) )
zpseu_nc = zpseur_nc
else
allocate( zpseu(2,norb,norb) )
zpseu = zpseur
endif
endif
!
! To find indices of initial and final slab
!
do ik=1, nrz
if (z(ik).le.zin+eps) kin=ik
if (z(ik).le.zfin-eps) kfin=ik
enddo
dz=z(kin+1)-z(kin)
allocate( amat( (2*n2d+npol*norb), (2*n2d+npol*norb) ) )
allocate( bmat( (2*n2d+npol*norb), (2*n2d+npol*norb) ) )
allocate( vec( (2*n2d+npol*norb), 2*(n2d+npol*nocros) ) )
allocate( aux( n2d, 2*n2d+npol*norb))
if (noncolin) then
allocate( zps_nc( norb*npol, norb*npol ) )
else
allocate( zps( norb, norb ) )
endif
amat=(0.d0,0.d0)
bmat=(0.d0,0.d0)
!
! zps=zpseu-e*qq for US-PP and zps=zpseu for norm-conserv. PP
!
orbin=orbin-1
do iorb=1, norbnow
do iorb1=1, norbnow
do iorb=1, norb
do iorb1=1, norb
if (noncolin) then
ij=0
do is=1,npol
do js=1,npol
ij=ij+1
zps_nc(npol*(iorb-1)+is, npol*(iorb1-1)+js)= &
zpseu_nc(orbin+iorb, orbin+iorb1,ij)
zps_nc(npol*(iorb-1)+is, npol*(iorb1-1)+js)= &
zpseu_nc(1,iorb,iorb1,ij)
if (lspinorb) then
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*zpseu_nc(2,iorb,iorb1,ij)
else
if (is.eq.js) &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*zpseu_nc(2,iorb,iorb1,ij)
endif
enddo
enddo
else
zps(iorb, iorb1)=zpseu(orbin+iorb,orbin+iorb1,iofspin)
zps(iorb,iorb1)=zpseu(1,iorb,iorb1)-eryd*zpseu(2,iorb,iorb1)
endif
enddo
enddo
do iorb=1, nocros+noins
nt=itnew(orbin+iorb)
if(tvanp(nt).or.lspinorb) then
ih=natih(orbin+iorb,2)
do iorb1=1, nocros+noins
if (natih(orbin+iorb,1).eq.natih(orbin+iorb1,1)) then
ih1=natih(orbin+iorb1,2)
if (noncolin) then
ij=0
do is=1,npol
do js=1,npol
ij=ij+1
if (lspinorb) then
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*qq_so(ih,ih1,ij,nt)
else
if (is.eq.js) &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*qq(ih,ih1,nt)
endif
enddo
enddo
else
zps(iorb,iorb1)=zps(iorb,iorb1)-eryd*qq(ih,ih1,nt)
endif
endif
enddo
endif
enddo
do iorb=1, nocros*npol
do iorb1=1, nocros*npol
if (noncolin) then
zps_nc(iorb+npol*(noins+nocros), iorb1+npol*(noins+nocros))= &
zps_nc(iorb,iorb1)
else
zps(iorb+noins+nocros,iorb1+noins+nocros)=zps(iorb,iorb1)
endif
enddo
enddo
!
! Forming the matrices A and B for generalized eigenvalue problem
@ -144,10 +125,11 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
bmat(ig+n2d,n)=fund0(ig, n)
enddo
enddo
!
! 2
!
do iorb=1, norbnow*npol
do iorb=1, norb*npol
do ig=1, n2d
amat(ig, 2*n2d+iorb)=funl1(ig, iorb)
amat(n2d+ig, 2*n2d+iorb)=fundl1(ig, iorb)
@ -158,13 +140,13 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
!
! 3
!
do iorb=1, norbnow*npol
do iorb=1, norb*npol
aorb=iorb
borb=iorb
if (iorb.le.npol*nocros) aorb=iorb+npol*(noins+nocros)
if (iorb.gt.npol*nocros) borb=iorb-npol*(noins+nocros)
do n=1, 2*n2d
do iorb1=1, norbnow*npol
do iorb1=1, norb*npol
if (noncolin) then
amat(2*n2d+iorb,n)=amat(2*n2d+iorb,n)+ &
zps_nc(aorb, iorb1)*intw1(iorb1,n)
@ -183,8 +165,8 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
! 4
!
do iorb=1, nocros*npol
do iorb1=1, norbnow*npol
do iorb2=1, norbnow*npol
do iorb1=1, norb*npol
do iorb2=1, norb*npol
if (noncolin) then
bmat(2*n2d+iorb,2*n2d+iorb1)=bmat(2*n2d+iorb,2*n2d+iorb1) &
-zps_nc(iorb,iorb2)*intw2(iorb2, iorb1)
@ -202,11 +184,11 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
!
! 5
!
do iorb=1, norbnow*npol
do iorb=1, norb*npol
aorb=iorb
if (iorb.le.npol*nocros) aorb=iorb+npol*(noins+nocros)
do iorb1=1, norbnow*npol
do iorb2=1, norbnow*npol
do iorb1=1, norb*npol
do iorb2=1, norb*npol
if (noncolin) then
amat(2*n2d+iorb,2*n2d+iorb1)=amat(2*n2d+iorb,2*n2d+iorb1)+ &
zps_nc(aorb,iorb2)*intw2(iorb2, iorb1)
@ -223,40 +205,27 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
!
! To reduce matrices and solve GEP A X = c B X; X = {a_n, a_\alpha}
!
call compbs_2(npol*nocros, npol*norbnow, n2d, 2*(n2d+npol*nocros), &
call compbs_2(npol*nocros, npol*norb, n2d, 2*(n2d+npol*nocros), &
amat, bmat, vec, kval, llapack)
!
! To normalize (over XY plane) all the states
!
! kfun=(0.d0,0.d0)
! kfund=(0.d0,0.d0)
! do ik=1, 2*(n2d+npol*nocros)
! do ig=1, n2d
! do j=1, 2*n2d+npol*norbnow
! kfun(ig,ik)= kfun(ig,ik)+amat(ig,j)*vec(j,ik)
! kfund(ig,ik)= kfund(ig,ik)+amat(n2d+ig,j)*vec(j,ik)
! enddo
! enddo
! enddo
!
call ZGEMM('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norbnow, &
one, amat, 2*n2d+npol*norbnow, vec, 2*n2d+npol*norbnow, &
call ZGEMM('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb, &
one, amat, 2*n2d+npol*norb, vec, 2*n2d+npol*norb, &
zero, kfun, n2d)
do ig=1,n2d
do ik=1, 2*n2d+npol*norbnow
do ik=1, 2*n2d+npol*norb
aux(ig,ik)=amat(n2d+ig,ik)
enddo
enddo
call ZGEMM('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norbnow, &
one, aux, n2d, vec, 2*n2d+npol*norbnow, zero, kfund, n2d)
call ZGEMM('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb, &
one, aux, n2d, vec, 2*n2d+npol*norb, zero, kfund, n2d)
do ik=1, 2*(n2d+npol*nocros)
raux=DDOT(2*n2d,kfun(1,ik),1,kfun(1,ik),1)*sarea
raux=1.d0/sqrt(raux)
call DSCAL(2*(2*n2d+npol*norbnow),raux,vec(1,ik),1)
call DSCAL(2*(2*n2d+npol*norb),raux,vec(1,ik),1)
call DSCAL(2*n2d,raux,kfun(1,ik),1)
call DSCAL(2*n2d,raux,kfund(1,ik),1)
enddo
@ -267,8 +236,8 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
call kbloch (2*(n2d+npol*nocros), kval)
call jbloch(2*(n2d+npol*nocros), n2d, norbf, norbnow, nocros, &
kfun, kfund, vec, kval, intw1, intw2, sarea, nchan, npol)
call jbloch(2*(n2d+npol*nocros), n2d, norbf, norb, nocros, &
kfun, kfund, vec, kval, intw1, intw2, nchan, npol)
!
! To save band structure result
!
@ -278,14 +247,14 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
!
! To account for the case of the right lead
!
if (lright.eq.1) then
if (lleft.eq.0) then
do i=1, 2*n2d
do j=1, 2*n2d+npol*norbnow
do j=1, 2*n2d+npol*norb
amat(i,j)=bmat(i,j)
enddo
enddo
do i=2*n2d+1, 2*n2d+npol*nocros
do j=1, 2*n2d+npol*norbnow
do j=1, 2*n2d+npol*norb
amat(i,j)=-bmat(i+npol*(nocros+noins),j)
enddo
enddo
@ -293,34 +262,16 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
!
! psi_k and psi'_k on the scattering region boundary
!
! do ik=1, 2*(n2d+npol*nocros)
! do ig=1, n2d
! do j=1, 2*n2d+npol*norbnow
! kfun(ig,ik)= kfun(ig,ik)+amat(ig,j)*vec(j,ik)
! kfund(ig,ik)= kfund(ig,ik)+amat(n2d+ig,j)*vec(j,ik)
! enddo
! enddo
! enddo
call ZGEMM('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norbnow, &
one, amat, 2*n2d+npol*norbnow, vec, 2*n2d+npol*norbnow, &
call ZGEMM('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb,&
one, amat, 2*n2d+npol*norb, vec, 2*n2d+npol*norb, &
zero, kfun, n2d)
do ig=1,n2d
do ik=1, 2*n2d+npol*norbnow
do ik=1, 2*n2d+npol*norb
aux(ig,ik)=amat(n2d+ig,ik)
enddo
enddo
call ZGEMM('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norbnow, &
one, aux, n2d, vec, 2*n2d+npol*norbnow, zero, kfund, n2d)
! do j=1,2*nchan
! write(6,'("------------------------------",i5,2f15.6)') j, kval(j)
! do ik=1,n2d
! if (ik.eq.n2d/2+1.and.noncolin) write(6,'("-------")')
! write(6,'(i5,f15.7)') ik, abs(kfun(ik,j))
! enddo
! enddo
! stop
call ZGEMM('n', 'n', n2d, 2*(n2d+npol*nocros), 2*n2d+npol*norb,&
one, aux, n2d, vec, 2*n2d+npol*norb, zero, kfund, n2d)
!
! kint(iorb, ik)=\sum_{iorb1} D_{iorb,iorb1}
@ -329,7 +280,7 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
!
do ik=1, 2*(n2d+npol*nocros)
do iorb=1, nocros*npol
do j=1, 2*n2d+npol*norbnow
do j=1, 2*n2d+npol*norb
kint(iorb,ik)=kint(iorb,ik)+amat(2*n2d+iorb,j)*vec(j,ik)
enddo
enddo
@ -341,7 +292,7 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
!
do ik=1, 2*(n2d+npol*nocros)
do iorb=1, nocros*npol
if (lright.eq.1) then
if (lleft.eq.0) then
kcoef(iorb,ik)=vec(2*n2d+iorb,ik)
else
kcoef(iorb,ik)=vec(2*n2d+npol*(nocros+noins)+iorb,ik)
@ -352,7 +303,7 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
!
! to set up B.S. for the right lead in the case of identical tips
!
if(lright.eq.0.and.ikind.eq.1) then
if(lleft.eq.1.and.ikind.eq.1) then
nchanr=nchan
call DCOPY(2*(n2d+npol*nocros), kval, 1, kvalr, 1)
kfunr=(0.d0,0.d0)
@ -360,19 +311,19 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
kintr=(0.d0,0.d0)
do i=1, 2*n2d
do j=1, 2*n2d+npol*norbnow
do j=1, 2*n2d+npol*norb
amat(i,j)=bmat(i,j)
enddo
enddo
do i=2*n2d+1, 2*n2d+npol*nocros
do j=1, 2*n2d+npol*norbnow
do j=1, 2*n2d+npol*norb
amat(i,j)=-bmat(i+npol*(nocros+noins),j)
enddo
enddo
do ik=1, 2*(n2d+npol*nocros)
do ig=1, n2d
do j=1, 2*n2d+npol*norbnow
do j=1, 2*n2d+npol*norb
kfunr(ig,ik)= kfunr(ig,ik)+amat(ig,j)*vec(j,ik)
kfundr(ig,ik)= kfundr(ig,ik)+amat(n2d+ig,j)*vec(j,ik)
enddo
@ -380,7 +331,7 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
enddo
do ik=1, 2*(n2d+npol*nocros)
do iorb=1, nocros*npol
do j=1, 2*n2d+npol*norbnow
do j=1, 2*n2d+npol*norb
kintr(iorb,ik)=kintr(iorb,ik)+amat(2*n2d+iorb,j)*vec(j,ik)
enddo
enddo
@ -396,9 +347,14 @@ subroutine compbs(lright, zin, zfin, nocros, norbnow, orbin, &
deallocate(amat)
deallocate(bmat)
deallocate(vec)
deallocate(zps)
deallocate(aux)
if (noncolin) deallocate(zps_nc)
if (noncolin) then
deallocate(zpseu_nc)
deallocate(zps_nc)
else
deallocate(zpseu)
deallocate(zps)
endif
call stop_clock('compbs')
return

23
PWCOND/compbs_2.f90
View File

@ -5,8 +5,9 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
subroutine compbs_2(nocros, norbnow, n2d, ntot, amat, bmat, &
subroutine compbs_2(nocros, norb, n2d, ntot, amat, bmat, &
vec, kval, llapack)
!
! This subroutine reduces amat, bmat into amt, bmt
! excluding inside lying orbitals and solves GEP:
@ -20,20 +21,20 @@ subroutine compbs_2(nocros, norbnow, n2d, ntot, amat, bmat, &
integer :: n2d, & ! 2D dimension
noins, & ! interior orbitals
nocros, & ! crossing orbitals
norbnow, & ! total number of orbitals
norb, & ! total number of orbitals
ntot ! ntot = 2*(n2d+nocros)
integer :: info, ishift, i, j, k, l
integer, allocatable :: ipiv(:)
logical :: llapack
complex(kind=DP) :: amat(2*n2d+norbnow, 2*n2d+norbnow), &
bmat(2*n2d+norbnow, 2*n2d+norbnow), &
vec(2*n2d+norbnow, ntot), kval(ntot)
complex(kind=DP) :: amat(2*n2d+norb, 2*n2d+norb), &
bmat(2*n2d+norb, 2*n2d+norb), &
vec(2*n2d+norb, ntot), kval(ntot)
complex(kind=DP), allocatable :: amt(:,:), bmt(:,:), &
auxa(:), auxb(:), auxc(:), &
hmat(:,:), hmt(:,:), vecaux(:,:)
call start_clock('compbs_2')
noins = norbnow-2*nocros
noins = norb-2*nocros
allocate( bmt( ntot, ntot ) )
allocate( amt( ntot, ntot ) )
@ -49,7 +50,7 @@ subroutine compbs_2(nocros, norbnow, n2d, ntot, amat, bmat, &
! To interchange inside and right crossing orbitals in a and b
!
! rows
do j=1, 2*n2d+norbnow
do j=1, 2*n2d+norb
do i=1, nocros
auxa(i)=amat(2*n2d+nocros+noins+i,j)
auxb(i)=bmat(2*n2d+nocros+noins+i,j)
@ -65,7 +66,7 @@ subroutine compbs_2(nocros, norbnow, n2d, ntot, amat, bmat, &
enddo
enddo
! columns
do j=1, 2*n2d+norbnow
do j=1, 2*n2d+norb
do i=1, nocros
auxa(i)=amat(j,2*n2d+nocros+noins+i)
auxb(i)=bmat(j,2*n2d+nocros+noins+i)
@ -135,7 +136,7 @@ subroutine compbs_2(nocros, norbnow, n2d, ntot, amat, bmat, &
!ccccccc
!
! Forming (2*n2d+norbnow, ntot) matrix of eigenvectors
! Forming (2*n2d+norb, ntot) matrix of eigenvectors
! cooficients, storing them in vec
!
vec=(0.d0,0.d0)
@ -165,7 +166,7 @@ subroutine compbs_2(nocros, norbnow, n2d, ntot, amat, bmat, &
! (to have a right order of orbitals again)
! rows
do j=1, 2*n2d+norbnow
do j=1, 2*n2d+norb
do i=1, nocros
auxa(i)=amat(2*n2d+nocros+i,j)
auxb(i)=bmat(2*n2d+nocros+i,j)
@ -181,7 +182,7 @@ subroutine compbs_2(nocros, norbnow, n2d, ntot, amat, bmat, &
enddo
enddo
! columns
do j=1, 2*n2d+norbnow
do j=1, 2*n2d+norb
do i=1, nocros
auxa(i)=amat(j,2*n2d+nocros+i)
auxb(i)=bmat(j,2*n2d+nocros+i)

136
PWCOND/cond_out.f90 Normal file
View File

@ -0,0 +1,136 @@
!
! 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 cond_out ()
use io_global, only : stdout
USE ions_base, ONLY: atm
use lsda_mod, only: nspin
USE noncollin_module, ONLY : noncolin, npol
use spin_orb, only : lspinorb
use cond
!---------------------------
! Some output
!---------------------------
implicit none
integer :: iorb, ipol, k
write(stdout,'(''----- General information -----'')')
write(stdout,*)
if(ikind.eq.0) then
write(stdout,'(''----- Complex band structure calculation -----'')')
elseif(ikind.eq.1) then
write(stdout,'(''--- T calc. with identical leads (ikind=1) --- '')')
elseif(ikind.eq.2) then
write(stdout,'(''--- T calc. with different leads (ikind=2) --- '')')
endif
if(nspin.eq.2) then
write(6,'(/,9x, ''LSDA calculations, spin index ='',i6)') iofspin
endif
if(nspin.eq.4) then
write(6,'(/,9x, ''Noncollinear calculations'')')
if(lspinorb) &
write(6,'(/,9x, ''Noncollinear calculations with spin-orbit'')')
endif
write (6, 300) nrx, nry, nz1
300 format (/,5x, &
& 'nrx = ',i12,/,5x, &
& 'nry = ',i12,/,5x, &
& 'nz1 = ',i12,/,5x)
write (6, 301) energy0, denergy, nenergy, ecut2d, ewind, epsproj
301 format (/,5x, &
& 'energy0 = ',1pe15.1,/,5x, &
& 'denergy = ',1pe15.1,/,5x, &
& 'nenergy = ',i10,/,5x, &
& 'ecut2d = ',1pe15.1,/,5x, &
& 'ewind = ',1pe15.1,/,5x, &
& 'epsproj = ',1pe15.1,/,5x)
if(ikind.eq.1) then
write(stdout,'(''----- Information about left/right lead -----'')')
else
write(stdout,'(''----- Information about left lead ----- '')')
endif
write (6, 200) nocrosl, noinsl, norbl, norbf, nrzl
200 format (/,5x, &
& 'nocros = ',i12,/,5x, &
& 'noins = ',i12,/,5x, &
& 'norb = ',i12,/,5x, &
& 'norbf = ',i12,/,5x, &
& 'nrz = ',i12,/,5x)
write(6, '(6x,''iorb type ibeta ang. mom.'',3x, &
& ''m position (a_0)'')')
write(6,'(5x,i4,4x,i5,5x,i3,6x,i3,6x,i3,'' taunew('', &
& i4,'')=('',3f8.4,'')'')') &
& ( iorb,tblml(1,iorb), tblml(2,iorb), tblml(3,iorb),&
& tblml(4,iorb), iorb, &
& (taunewl(ipol,iorb),ipol=1,3), iorb=1, norbl )
if(norbl.le.80) then
write(6,'(4x,''k slab'',3x,'' z(k) z(k+1)'', &
& 5x,''crossing(iorb=1,norb)'')')
do k=1, nrzl
write(6,'(2x,i3,2x,3f7.4,3x,80i1)') &
k,zl(k),zl(k+1),zl(k+1)-zl(k),(crosl(iorb,k),iorb=1,norbl)
enddo
endif
if(ikind.eq.2) then
write(stdout,'(''----- Information about right lead -----'')')
write (6, 200) nocrosr, noinsr, norbr, norbf, nrzr
write(6, '(6x,''iorb type ibeta ang. mom.'',3x, &
& ''m position (a_0)'')')
write(6,'(5x,i4,4x,i5,5x,i3,6x,i3,6x,i3,'' taunew('', &
& i4,'')=('',3f8.4,'')'')') &
& ( iorb,tblmr(1,iorb), tblmr(2,iorb), tblmr(3,iorb),&
& tblmr(4,iorb), iorb, &
& (taunewr(ipol,iorb),ipol=1,3), iorb=1, norbr )
if(norbr.le.80) then
write(6,'(4x,''k slab'',3x,'' z(k) z(k+1)'', &
& 5x,''crossing(iorb=1,norb)'')')
do k=1, nrzr
write(6,'(2x,i3,2x,3f7.4,3x,80i1)') &
k,zr(k),zr(k+1),zr(k+1)-zr(k),(crosr(iorb,k),iorb=1,norbr)
enddo
endif
endif
if(ikind.gt.0) then
write(6,'(''----- Information about scattering region -----'')')
write (6, 201) noinss, norbs, norbf, nrzs
201 format (/,5x, &
& 'noins = ',i12,/,5x, &
& 'norb = ',i12,/,5x, &
& 'norbf = ',i12,/,5x, &
& 'nrz = ',i12,/,5x)
write(6, '(6x,''iorb type ibeta ang. mom.'',3x, &
& ''m position (a_0)'')')
write(6,'(5x,i4,4x,i5,5x,i3,6x,i3,6x,i3,'' taunew('', &
& i4,'')=('',3f8.4,'')'')') &
& ( iorb,tblms(1,iorb), tblms(2,iorb), tblms(3,iorb),&
& tblms(4,iorb), iorb, &
& (taunews(ipol,iorb),ipol=1,3), iorb=1, norbs )
if(norbs.le.80) then
write(6,'(4x,''k slab'',3x,'' z(k) z(k+1)'', &
& 5x,''crossing(iorb=1,norb)'')')
do k=1, nrzs
write(6,'(2x,i3,2x,3f7.4,3x,80i1)') &
k,zs(k),zs(k+1),zs(k+1)-zs(k),(cross(iorb,k),iorb=1,norbs)
enddo
endif
endif
return
end subroutine cond_out

253
PWCOND/condcom.f90
View File

@ -1,4 +1,3 @@
!
! Copyright (C) 2003 A. Smogunov
! This file is distributed under the terms of the
@ -12,139 +11,150 @@
!
!
!
MODULE zdir_cond
USE kinds
!
! ... description of z direction
MODULE geomcell_cond
USE kinds, only : DP
!
SAVE
!
INTEGER :: &
nrz, &! total number of the slabs in z direction
nrzp, &! number of slabs per CPU
dnslab, &! additional slabs (nrz=nr3s+dnslab)
nz1 ! number of subslabs in the slab
INTEGER :: &
nrx, & ! number of mesh points in the x direction
nry, & ! -||- y direction
nrzl, & ! number of slabsh in the z direction for the left lead
nrzs, & ! -||- for the scatt. region
nrzr, & ! -||- for the right lead
nrzpl, & ! number of slabs per CPU for the left lead
nrzps, & ! -||- for the scatt. region
nrzpr, & ! -||- for the right lead
ngper, & ! number of perpendicular G vectors
ngpsh, & ! number of shells for G
nkpts, & ! number of kpts in the perpendicular direction
n2d, & ! dimension of reduced vector space in XY
nz1 ! number of subslabs in the slab
INTEGER, ALLOCATABLE :: &
nkofz(:) ! local CPU order of the slabs
REAL(KIND=DP) :: &
bdl1, bdl2, &! boundaries for the left tip
bdr1, bdr2, &! --- || --- right tip
zl ! the lenght of the unit cell in z direction
REAL(KIND=DP), ALLOCATABLE :: &
z(:) ! the unit cell division in z direction
!
END MODULE zdir_cond
ninsh(:) ! number of G in shell
REAL(kind=DP) :: &
bdl, & ! right boundary of the left lead
bds, & ! -||- of the scatt. region
bdr, & ! -||- of the right lead
sarea ! the cross section
REAL(kind=DP), ALLOCATABLE :: &
zl(:), & ! the division in the z direction of the left lead
zs(:), & ! -||- of the scatt. reg.
zr(:), & ! -||- of the right lead
xyk(:,:), & ! coordinates of perpendicular k points
wkpt(:), & ! the weight of k point
gper(:,:),& ! coordinates of perpendicular G
gnsh(:) ! the norm of the G shell
END MODULE geomcell_cond
!
!
MODULE perp_cond
USE kinds
!
! ... description of the direction perpendicular to z (XY)
MODULE orbcell_cond
USE parameters, only : ndmx, nbrx, npsx
USE kinds, only : DP
!
! description of nonlocal orbitals
SAVE
!
INTEGER :: &
nrx, &! the number of mesh points in x direction
nry, &! -- in y direction
ngper, &! number of perpendicular G vectors
ngpsh, &! number of shells for G
nkpts, &! number of kpts in the perpendicular direction
n2d ! dimension of reduced vector space in XY
INTEGER :: &
norbl, & ! number of orbitals for the left lead
norbs, & ! -||- for the scatt. region
norbr, & ! -||- for the right lead
nocrosl, & ! number of crossing orbitals for left lead
nocrosr, & ! -||- for the right lead
noinsl, & ! number of interior orbitals for the left lead
noinss, & ! -||- for the scatt. region
noinsr, & ! -||- for the right lead
norbf ! max number of needed orbitals
INTEGER, ALLOCATABLE :: &
ninsh(:) ! number of G in shell
REAL(KIND=DP) :: &
sarea ! the cross section
REAL(KIND=DP), ALLOCATABLE :: &
xyk(:,:), &! coordinates of perpendicular k point
wkpt(:), &
gper(:,:), &! coordinates of perpendicular G
gnsh(:) ! the norm of the shell
!
END MODULE perp_cond
tblml(:,:), & ! the type/beta/l/m of each orbital for the left lead
tblms(:,:), & ! -||- for the scatt. reg.
tblmr(:,:), & ! -||- for the right lead
crosl(:,:), & ! 1 if the orbital crosses the slab - for the left lead
cross(:,:), & ! -||- for the scatt. reg.
crosr(:,:) ! -||- for the right lead
REAL(kind=DP) :: &
rl(ndmx,npsx), & ! radial mesh for the left lead
rs(ndmx,npsx), & ! -||- for the scatt. reg.
rr(ndmx,npsx), & ! -||- for the right lead
rabl(ndmx,npsx), & ! log. mesh for the left lead
rabs(ndmx,npsx), & ! -||- for the scatt. reg.
rabr(ndmx,npsx), & ! -||- for the right lead
betarl(ndmx,nbrx,npsx), & ! beta functions for the left lead
betars(ndmx,nbrx,npsx), & ! -||- for the scatt. reg.
betarr(ndmx,nbrx,npsx) ! -||- for the right lead
REAL(kind=DP), ALLOCATABLE :: &
taunewl(:,:), & ! center of each orbital and its radius - left lead
taunews(:,:), & ! -||- - scatt. reg.
taunewr(:,:), & ! -||- - right lead
zpseul(:,:,:), & ! coefficients of nonlocal pseudopotential - left lead
zpseus(:,:,:), & ! -||- - scatt. reg.
zpseur(:,:,:) ! -||- - right lead
COMPLEX(kind=DP), ALLOCATABLE :: &
zpseul_nc(:,:,:,:), &! coefficients of nonlocal PP (nc case) - left lead
zpseus_nc(:,:,:,:), &! -||- - scatt. reg.
zpseur_nc(:,:,:,:) ! -||- - right lead
END MODULE orbcell_cond
!
!
MODULE eigen_cond
USE kinds
!
! ... Eigenvalue equation for local potential
USE kinds, only : DP
!
! Eigenvalue equation for local potential
SAVE
!
COMPLEX(KIND=DP), ALLOCATABLE :: &
vppot(:,:,:,:), &! Fourier comp. of local potential in each slab
psiper(:,:,:), &! eigenvectors in each slab
newbg(:,:), &! basis of reduced set --> exp(G)
zk(:,:) ! the k for each eigenvalue (computed through zkr)
REAL(KIND=DP), ALLOCATABLE :: &
zkr(:,:) ! eigenvalues for first energy point
!
END MODULE eigen_cond
COMPLEX(kind=DP), ALLOCATABLE :: &
vppotl(:,:,:,:), & ! Fourier comp. of local potential in each slab - left lead
vppots(:,:,:,:), & ! -||- - scatt. reg.
vppotr(:,:,:,:), & ! -||- - right lead
psiperl(:,:,:), & ! eigenvectors in each slab - left lead
psipers(:,:,:), & ! -||- - scatt. reg.
psiperr(:,:,:), & ! -||- - right lead
zkl(:,:), & ! the k for each eigenvalue (computed through zkr) - left lead
zks(:,:), & ! -||- - scatt. reg.
zkr(:,:), & ! -||- - right lead
newbg(:,:) ! reduced basis set --> exp(G)
REAL(kind=DP), ALLOCATABLE :: &
zkrl(:,:), & ! 2d eigenvalues - left lead
zkrs(:,:), & ! -||- - scatt. reg.
zkrr(:,:) ! -||- - right lead
END MODULE eigen_cond
!
!
MODULE control_cond
USE kinds
!
! ... control of the run
USE kinds, only : DP
!
! control of the run
SAVE
!
INTEGER :: &
ikind, &! the kind of calculation
nenergy, &! number of energies computed
iofspin ! spin index for calculation
REAL(KIND=DP) :: &
energy0, &! initial energy
denergy, &! delta of energy
eryd, &! current energy
ecut2d, &! 2D cutoff
ewind, &! the window above energy for 2D computation
delgep, &! infinitesimal for GEP
epsproj, &! accuracy of n2d reduction
cutplot ! cutoff of Im(k) for CB plotting
REAL(KIND=DP), ALLOCATABLE :: &
earr(:), &! energy array
tran_tot(:) ! transmission array
LOGICAL :: &
lwrite_loc, &! if .T. save eigenproblem result on fil_loc
lread_loc, &! if .T. read eigenproblem result from fil_loc
llapack, &! if .T. use LAPACK routine for GEP
llocal ! if .T. the local implementation
!
INTEGER :: &
orbj_in, orbj_fin, &
ikind, & ! the kind of calculation
nenergy, & ! number of energies computed
iofspin ! spin index for calculation
REAL(kind=DP) :: &
efl, & ! the Ef of the left lead
efs, & ! the Ef of the scatt. reg.
efr, & ! the Ef of the right lead
energy0, & ! initial energy
eryd, & ! the current energy in Ry
denergy, & ! delta of energy
ecut2d, & ! 2D cutoff
ewind, & ! the window above energy for 2D computation
delgep, & ! infinitesimal for GEP
epsproj, & ! accuracy of n2d reduction
cutplot ! cutoff of Im(k) for CB plotting
REAL(kind=DP), ALLOCATABLE :: &
earr(:), & ! energy array
tran_tot(:) ! transmission array
LOGICAL :: &
lwrite_loc, & ! if .t. save eigenproblem result on fil_loc
lread_loc, & ! if .t. read eigenproblem result from fil_loc
llapack, & ! if .t. use LAPACK routine for GEP
lwrite_cond, & ! if .t. save variables needed for pwcond
lread_cond, & ! if .t. read variables needed for pwcond
llocal ! if .t. the local implementation
END MODULE control_cond
!
!
MODULE cross_cond
USE kinds
!
! ... description of nonlocal orbitals
!
SAVE
!
INTEGER :: &
nocrosl, &! number of crossing orbitals for left tip
nocrosr, &! -- || -- for right tip
noinsl, &! number of interior orbitals for left tip
noinss, &! -- || -- for scat. reg.
noinsr, &! -- || -- for right tip
norb, &! total number of orbitals
norbf, &! max number of needed orbitals
norbs ! tot. number of orbitals in scat. region
INTEGER, ALLOCATABLE :: &
itnew(:), &! the type of each orbital
nbnew(:), &! the nb of each orbital
ls(:), &! the l of each orbital
mnew(:), &! the m of each orbital
cross(:,:), &! 1 if the orbital crosses the slab
natih(:,:) ! na and ih of the orbital
REAL(KIND=DP), ALLOCATABLE :: &
rsph(:,:), &! the radius of nonlocal sphere
taunew(:,:), &! center of each orbital
zpseu(:,:,:) ! coefficients of nonlocal pseudopotential
COMPLEX(KIND=DP), ALLOCATABLE :: &
zpseu_nc(:,:,:) ! coefficients of nonlocal pseudopotential (nc case)
!
END MODULE cross_cond
!
!
MODULE scattnl_cond
USE kinds
@ -171,22 +181,22 @@ END MODULE scattnl_cond
MODULE cb_cond
USE kinds
!
! ... Some variables of CBS for the tips needed for matching
! ... Some variables of CBS for the leads needed for matching
!
SAVE
!
INTEGER :: &
nchanl, &! number of prop. channels in the left tip
nchanr ! -- || -- right tip
nchanl, &! number of prop. channels in the left lead
nchanr ! -- || -- right lead
COMPLEX(KIND=DP), ALLOCATABLE :: &
kvall(:), &! k for left tip
kfunl(:,:), &! phi_k(z=d) for left tip
kfundl(:,:), &! phi_k'(z=d) for left tip
kvall(:), &! k for the left lead
kfunl(:,:), &! phi_k(z=d) for the left lead
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.
kvalr(:), &! k for right tip
kfunr(:,:), &! phi_k(z=0) for right tip
kfundr(:,:), &! phi_k'(z=0) for right tip
kvalr(:), &! k for the right lead
kfunr(:,:), &! phi_k(z=0) for the right lead
kfundr(:,:), &! phi_k'(z=0) for the right lead
kintr(:,:), &! integral of phi_k with beta-fun.
kcoefr(:,:) ! coeff. of phi_k over nonloc. fun.
!
@ -194,11 +204,10 @@ END MODULE cb_cond
!
!
MODULE cond
USE zdir_cond
USE perp_cond
use geomcell_cond
USE orbcell_cond
USE eigen_cond
USE control_cond
USE cross_cond
USE scattnl_cond
USE cb_cond
END MODULE cond

224
PWCOND/do_cond.f90
View File

@ -24,16 +24,15 @@ SUBROUTINE do_cond(nodenumber)
IMPLICIT NONE
CHARACTER(len=3) nodenumber
REAL(kind=DP), PARAMETER :: eps=1.d-8
REAL(kind=DP) :: wtot
INTEGER :: ik, ien, ios, orbin, orbfin
LOGICAL :: write0
INTEGER :: ik, ien, ios
NAMELIST /inputcond/ outdir, prefix, band_file, tran_file, fil_loc, &
lwrite_loc, lread_loc, ikind, iofspin, llocal, &
bdl1, bdl2, bdr1, bdr2, nz1, dnslab, energy0, &
denergy, ecut2d, ewind, epsproj, delgep,cutplot,&
llapack
NAMELIST /inputcond/ outdir, prefixt, prefixl, prefixs, prefixr, &
band_file, tran_file, save_file, fil_loc, &
lwrite_loc, lread_loc, lwrite_cond, lread_cond, &
orbj_in,orbj_fin,ikind,iofspin,llocal,llapack, &
bdl, bds, bdr, nz1, energy0, denergy, ecut2d, &
ewind, epsproj, delgep, cutplot
CHARACTER (LEN=80) :: input_file
INTEGER :: nargs, iiarg, ierr, ILEN
@ -47,28 +46,34 @@ SUBROUTINE do_cond(nodenumber)
! set default values for variables in namelist
!
outdir = './'
prefix = ' '
prefixt = ' '
prefixl = ' '
prefixs = ' '
prefixr = ' '
band_file = ' '
tran_file = ' '
save_file = ' '
fil_loc = ' '
lwrite_loc = .FALSE.
lread_loc = .FALSE.
lwrite_cond = .FALSE.
lread_cond = .FALSE.
orbj_in = 0
orbj_fin = 0
iofspin = 1
ikind = 0
bdl1 = 0.d0
bdl2 = 0.d0
bdr1 = 0.d0
bdr2 = 0.d0
bdl = 0.d0
bds = 0.d0
bdr = 0.d0
nz1 = 11
dnslab = 0
energy0 = 0.d0
denergy = 0.d0
ecut2d = 0.d0
ewind = 0.d0
ewind = 1.d0
llocal = .FALSE.
llapack = .FALSE.
llapack = .TRUE.
epsproj = 1.d-3
delgep = 0.d0
delgep = 5.d-10
cutplot = 2.d0
IF ( ionode ) THEN
@ -140,21 +145,27 @@ SUBROUTINE do_cond(nodenumber)
! ... Broadcast variables
!
CALL mp_bcast( tmp_dir, ionode_id )
CALL mp_bcast( prefix, ionode_id )
CALL mp_bcast( prefixt, ionode_id )
CALL mp_bcast( prefixl, ionode_id )
CALL mp_bcast( prefixs, ionode_id )
CALL mp_bcast( prefixr, ionode_id )
CALL mp_bcast( band_file, ionode_id )
CALL mp_bcast( tran_file, ionode_id )
CALL mp_bcast( fil_loc, ionode_id )
CALL mp_bcast( save_file, ionode_id )
CALL mp_bcast( lwrite_loc, ionode_id )
CALL mp_bcast( lread_loc, ionode_id )
CALL mp_bcast( lwrite_cond, ionode_id )
CALL mp_bcast( lread_cond, ionode_id )
CALL mp_bcast( ikind, ionode_id )
CALL mp_bcast( iofspin, ionode_id )
CALL mp_bcast( orbj_in, ionode_id )
CALL mp_bcast( orbj_fin, ionode_id )
CALL mp_bcast( llocal, ionode_id )
CALL mp_bcast( bdl1, ionode_id )
CALL mp_bcast( bdl2, ionode_id )
CALL mp_bcast( bdr1, ionode_id )
CALL mp_bcast( bdr2, ionode_id )
CALL mp_bcast( bdl, ionode_id )
CALL mp_bcast( bds, ionode_id )
CALL mp_bcast( bdr, ionode_id )
CALL mp_bcast( nz1, ionode_id )
CALL mp_bcast( dnslab, ionode_id )
CALL mp_bcast( energy0, ionode_id )
CALL mp_bcast( denergy, ionode_id )
CALL mp_bcast( ecut2d, ionode_id )
@ -177,75 +188,134 @@ SUBROUTINE do_cond(nodenumber)
CALL mp_bcast( earr, ionode_id )
CALL mp_bcast( tran_tot, ionode_id )
!
! Now allocate space for pwscf variables, read and check them.
!
CALL read_file
CALL openfil
CALL struc_fact (nat,tau,ntyp,ityp,ngm,g,bg, &
nr1,nr2,nr3,strf,eigts1,eigts2,eigts3)
CALL init_us_1
CALL newd
IF (lread_cond) THEN
call save_cond (.false.,1,efl,nrzl,nocrosl,noinsl, &
norbl,rl,rabl,betarl)
if(ikind.eq.1) then
call save_cond (.false.,2,efs,nrzs,-1, &
noinss,norbs,rs,rabs,betars)
norbr = norbl
nocrosr = nocrosl
noinsr = noinsl
endif
if(ikind.eq.2) then
call save_cond (.false.,2,efs,nrzs,-1, &
noinss,norbs,rs,rabs,betars)
call save_cond (.false.,3,efr,nrzr,nocrosr,&
noinsr,norbr,rr,rabr,betarr)
endif
ELSE
IF (prefixt.ne.' ') then
prefix = prefixt
call read_file
call init_us_1
call newd
IF (ikind.eq.0) then
CALL init_cond(1,'t')
ELSEIF (ikind.eq.1) then
CALL init_cond(2,'t')
ELSEIF (ikind.eq.2) then
CALL init_cond(3,'t')
ENDIF
CALL clean_pw(.true.)
ENDIF
IF (prefixl.ne.' ') then
prefix = prefixl
call read_file
call init_us_1
call newd
CALL init_cond(1,'l')
CALL clean_pw(.true.)
ENDIF
IF (prefixs.ne.' ') then
prefix = prefixs
call read_file
call init_us_1
call newd
CALL init_cond(1,'s')
CALL clean_pw(.true.)
ENDIF
IF (prefixr.ne.' ') then
prefix = prefixr
call read_file
call init_us_1
call newd
CALL init_cond(1,'r')
CALL clean_pw(.true.)
ENDIF
ENDIF
IF (lwrite_cond) then
call save_cond (.true.,1,efl,nrzl,nocrosl,noinsl, &
norbl,rl,rabl,betarl)
if(ikind.gt.0) call save_cond (.true.,2,efs,nrzs,-1, &
noinss,norbs,rs,rabs,betars)
if(ikind.gt.1) call save_cond (.true.,3,efr,nrzr,nocrosr,&
noinsr,norbr,rr,rabr,betarr)
endif
call cond_out
!
! Allocation for pwcond variables
!
CALL allocate_cond
CALL init_cond
CALL stop_clock('init')
IF (llocal) &
CALL local_set(nocrosl,noinsl,noinss,nocrosr,noinsr,norb,norbs)
CALL poten
CALL local_set(nocrosl,noinsl,norbl,noinss,norbs,nocrosr,noinsr,norbr)
DO ik=1, nkpts
do ik=1, nkpts
CALL init_gper(ik)
!
! The main loop
!
eryd = earr(1)/rytoev + ef
CALL local
DO ien=1, nenergy
eryd = earr(ien)/rytoev + ef
CALL form_zk(n2d, nrzp, zkr, zk, eryd, tpiba)
orbin=1
orbfin=orbin-1+2*nocrosl+noinsl
CALL scatter_forw(bdl1, bdl2, orbin, orbfin)
orbin=1
orbfin=orbin-1+2*nocrosl+noinsl
CALL compbs(0, bdl1, bdl2, nocrosl, &
orbfin-orbin+1, orbin, nchanl, kvall, &
kfunl, kfundl, kintl, kcoefl)
IF (ikind.EQ.2) THEN
orbin=2*nocrosl+noinsl+noinss+1
orbfin=orbin-1+2*nocrosr+noinsr
CALL scatter_forw(bdr1, bdr2, orbin, orbfin)
orbin=2*nocrosl+noinsl+noinss+1
orbfin=orbin-1+2*nocrosr+noinsr
CALL compbs(1, bdr1, bdr2, nocrosr, &
orbfin-orbin+1, orbin, nchanr, kvalr, &
kfunr, kfundr, kintr, kcoefr)
ENDIF
CALL summary_band(ik,ien)
IF (ikind.NE.0) THEN
orbin=nocrosl+noinsl+1
orbfin=orbin-1+nocrosl+noinss+nocrosr
CALL scatter_forw(bdl2, bdr1, orbin, orbfin)
CALL transmit(ik,ien)
ENDIF
ENDDO
CALL free_mem
ENDDO
call local
do ien=1, nenergy
eryd = earr(ien)/rytoev + efl
call form_zk(n2d, nrzpl, zkrl, zkl, eryd, tpiba)
call scatter_forw(nrzl, nrzpl, zl, psiperl, zkl, norbl, &
tblml, crosl, taunewl, rl, rabl, betarl)
call compbs(1, nocrosl, norbl, nchanl, kvall, kfunl, kfundl, &
kintl, kcoefl)
IF (ikind.EQ.2) THEN
eryd = earr(ien)/rytoev + efr
call form_zk(n2d, nrzpr, zkrr, zkr, eryd, tpiba)
call scatter_forw(nrzr, nrzpr, zr, psiperr, zkr, norbr, &
tblmr, crosr, taunewr, rr, rabr, betarr)
call compbs(0, nocrosr, norbr, nchanr, kvalr, kfunr, kfundr,&
kintr, kcoefr)
ENDIF
call summary_band(ik,ien)
IF (ikind.NE.0) THEN
eryd = earr(ien)/rytoev + efs
call form_zk(n2d, nrzps, zkrs, zks, eryd, tpiba)
call scatter_forw(nrzs, nrzps, zs, psipers, zks, norbs, &
tblms, cross, taunews, rs, rabs, betars)
write(6,*) 'to transmit'
call transmit(ik, ien)
endif
enddo
call free_mem
enddo
IF(ikind.GT.0.AND.tran_file.NE.' ') &
CALL summary_tran(tran_file,nenergy,earr,tran_tot)
CALL summary_tran()
CALL print_clock_pwcond()
CALL stop_clock('PWCOND')
RETURN
END SUBROUTINE do_cond
return
END SUBROUTINE do_cond

5
PWCOND/form_zk.f90
View File

@ -7,14 +7,13 @@
!
subroutine form_zk(n2d, nrzp, zkr, zk, e, tpiba)
!
! Just to construct complex zk=sqrt(e-E_n) for energy e
! from eigenvalues E_n found for some initial energy
! To construct complex wavevectors zk=sqrt(e-E_n)
! for an energy e from eigenvalues E_n of 2d problem
!
USE kinds, only : DP
implicit none
integer :: nrzp, n2d, n, k
real(kind=DP) :: zkr(n2d,nrzp), e, ed, tpiba
real(kind=DP), parameter :: eps=1.d-4
complex(kind=DP) :: zk(n2d,nrzp)
do k=1, nrzp

100
PWCOND/four.f90
View File

@ -1,4 +1,3 @@
!
! Copyright (C) 2003 A. Smogunov
! This file is distributed under the terms of the
@ -6,7 +5,7 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
subroutine four(alpha, w0, k, dz)
subroutine four(w0, z0, dz, tblm, taunew, r, rab, betar)
!
! This routine computes the bidimensional fourier transform of the
! beta function. It has been implemented for s, p, d-orbitals.
@ -22,29 +21,31 @@ subroutine four(alpha, w0, k, dz)
! p_z, p_{-x}, p_{-y} ;
! d_{z^2-1}, d_{-xz}, d_{-yz}, d_{x^2-y^2}, d_{xy}
!
! input: alpha - the number of the orbital in the full list.
! This orbital is the fist one for a given l.
! k - the number of the slab in z direction
! dz - the slab width
! input: tblm - array characterizing the orbital.
! taunew - coordinates and radius of the orbital.
! z0 - the initial z
! dz - the slab width
!
! output: w0(z, g, m), where
! z(k)< z <z(k+1) and
! z0< z <z0+dz
! g - 2D g-vector
!
#include "f_defs.h"
use pwcom
USE parameters, ONLY: ndmx
USE uspp_param, ONLY : betar
use atom, only: msh, r, rab
use cond
USE kinds, ONLY: DP
USE constants, ONLY : tpi, fpi
USE parameters, only : ndmx
USE cell_base, ONLY : alat, tpiba
USE cond, ONLY : sarea, nz1, ngper, gper, ninsh, gnsh, ngpsh
implicit none
integer :: alpha, k, kz, ig, ign, igphi, il, il1, nbb, itp, &
indexr, iz, lb, ir, nmesh, nmeshs
indexr, iz, lb, ir, nmesh, nmeshs, tblm(4)
real(kind=DP), parameter :: eps=1.d-8
complex(kind=DP), parameter :: cim=(0.d0, 1.d0)
real(kind=DP) :: gn, s1, s2, cs, sn, cs2, sn2, arg, rz, dz1, zr, &
dr, dz, bessj
dr, z0, dz, bessj, taunew(4), r(ndmx), &
rab(ndmx), betar(ndmx)
real(kind=DP), allocatable :: x1(:), x2(:), x3(:), x4(:)
real(kind=DP), allocatable :: fx1(:), fx2(:), fx3(:), fx4(:), zsl(:)
complex(kind=DP) :: exg, cong, xfact, w0(nz1, ngper, 5)
@ -62,89 +63,82 @@ implicit none
allocate( zsl( nz1) )
allocate( wadd( nz1, ngper ) )
itp=itnew(alpha)
nbb=nbnew(alpha)
lb=ls(alpha)
nmesh=indexr(rsph(nbb,itp)*alat,msh(itp),r(1,itp))
lb = tblm(3)
nmesh=indexr(taunew(4)*alat,ndmx,r)
dz1=dz/nz1
zsl(1)=(z(k)+dz1*0.5d0-taunew(3,alpha))*alat
do kz=2, nz1
zsl(kz)=zsl(kz-1)+dz1*alat
zsl(1)=(z0+dz1*0.5d0-taunew(3))*alat
do kz = 2, nz1
zsl(kz) = zsl(kz-1)+dz1*alat
enddo
ig=0
do ign=1, ngpsh
gn=gnsh(ign)
do kz=1, nz1
if (abs(zsl(kz))+eps.le.rsph(nbb,itp)*alat) then
iz=indexr(zsl(kz),nmesh,r(1,itp))
if (abs(zsl(kz))+eps.le.taunew(4)*alat) then
iz=indexr(zsl(kz),nmesh,r)
if ((nmesh-iz)/2*2.eq.nmesh-iz) then
nmeshs=nmesh
else
nmeshs=nmesh+1
endif
do ir=iz, nmeshs
rz=sqrt(r(ir,itp)**2-zsl(kz)**2)
rz=sqrt(r(ir)**2-zsl(kz)**2)
if (lb.eq.0) then
x1(ir)=betar(ir,nbb,itp)*bessj(0,gn*rz)
x1(ir)=betar(ir)*bessj(0,gn*rz)
elseif (lb.eq.1) then
x1(ir)=betar(ir,nbb,itp)*bessj(1,gn*rz)/ &
r(ir,itp)*rz
x2(ir)=betar(ir,nbb,itp)*bessj(0,gn*rz)/ &
r(ir,itp)
x1(ir)=betar(ir)*bessj(1,gn*rz)/r(ir)*rz
x2(ir)=betar(ir)*bessj(0,gn*rz)/r(ir)
elseif (lb.eq.2) then
x1(ir)=betar(ir,nbb,itp) &
*bessj(2,gn*rz)*rz**2/r(ir,itp)**2
x2(ir)=betar(ir,nbb,itp) &
*bessj(1,gn*rz)*rz/r(ir,itp)**2
x3(ir)=betar(ir,nbb,itp)*bessj(0,gn*rz)/ &
r(ir,itp)**2
x4(ir)=betar(ir,nbb,itp)*bessj(0,gn*rz)
x1(ir)=betar(ir)*bessj(2,gn*rz)*rz**2/r(ir)**2
x2(ir)=betar(ir)*bessj(1,gn*rz)*rz/r(ir)**2
x3(ir)=betar(ir)*bessj(0,gn*rz)/r(ir)**2
x4(ir)=betar(ir)*bessj(0,gn*rz)
else
call errore ('four','ls not programmed ',1)
endif
enddo
call simpson(nmeshs-iz+1,x1(iz),rab(iz,itp),fx1(kz))
call simpson(nmeshs-iz+1,x1(iz),rab(iz),fx1(kz))
if (iz.eq.1) then
dr=r(iz,itp)
dr=r(iz)
else
dr=r(iz,itp)-r(iz-1,itp)
dr=r(iz)-r(iz-1)
endif
zr=r(iz,itp)-abs(zsl(kz))
zr=r(iz)-abs(zsl(kz))
if (lb.eq.0) then
if (iz.eq.1) then
x1(iz-1)=betar(iz,nbb,itp)-betar(iz,nbb,itp)/dr*zr
x1(iz-1)=betar(iz)-betar(iz)/dr*zr
else
x1(iz-1)=betar(iz,nbb,itp)- &
(betar(iz,nbb,itp)-betar(iz-1,nbb,itp))/dr*zr
x1(iz-1)=betar(iz)-(betar(iz)-betar(iz-1))/dr*zr
endif
fx1(kz)=fx1(kz)+(x1(iz-1)+x1(iz))*0.5d0*zr
else
fx1(kz)=fx1(kz)+x1(iz)*0.5d0*zr
call simpson(nmeshs-iz+1,x2(iz),rab(iz,itp),fx2(kz))
call simpson(nmeshs-iz+1,x2(iz),rab(iz),fx2(kz))
endif
if (lb.eq.1) then
if(iz.eq.1) then
x2(iz-1)=0.d0
else
x2(iz-1)=(betar(iz,nbb,itp)-(betar(iz,nbb,itp)- &
betar(iz-1,nbb,itp))/dr*zr)/abs(zsl(kz))
x2(iz-1)=(betar(iz)-(betar(iz)- &
betar(iz-1))/dr*zr)/abs(zsl(kz))
endif
fx2(kz)=fx2(kz)+(x2(iz-1)+x2(iz))*0.5d0*zr
endif
if (lb.eq.2) then
fx2(kz)=fx2(kz)+x2(iz)*0.5d0*zr
call simpson(nmeshs-iz+1,x3(iz),rab(iz,itp),fx3(kz))
call simpson(nmeshs-iz+1,x4(iz),rab(iz,itp),fx4(kz))
call simpson(nmeshs-iz+1,x3(iz),rab(iz),fx3(kz))
call simpson(nmeshs-iz+1,x4(iz),rab(iz),fx4(kz))
if(iz.eq.1) then
x3(iz-1)=0.d0
x4(iz-1)=0.d0
else
x3(iz-1)=(betar(iz,nbb,itp)-(betar(iz,nbb,itp)- &
betar(iz-1,nbb,itp))/dr*zr)/abs(zsl(kz))**2
x4(iz-1)=betar(iz,nbb,itp)-(betar(iz,nbb,itp)- &
betar(iz-1,nbb,itp))/dr*zr
x3(iz-1)=(betar(iz)-(betar(iz)- &
betar(iz-1))/dr*zr)/abs(zsl(kz))**2
x4(iz-1)=betar(iz)-(betar(iz)- &
betar(iz-1))/dr*zr
endif
fx3(kz)=fx3(kz)+(x3(iz-1)+x3(iz))*0.5d0*zr
fx4(kz)=fx4(kz)+(x4(iz-1)+x4(iz))*0.5d0*zr
@ -237,7 +231,7 @@ function indexr(zz, ndim, r)
!
! abs(zz)<r(indexr)
!
iz=1
iz = 1
do while(r(iz).le.abs(zz)+1.d-10)
iz=iz+1
enddo

16
PWCOND/free_mem.f90
View File

@ -13,14 +13,18 @@ subroutine free_mem
use cond
implicit none
!
! From local_2
! From allocate_cond
!
deallocate(psiper)
deallocate(zk)
deallocate(psiperl)
deallocate(zkl)
deallocate(zkrl)
deallocate(psipers)
deallocate(zks)
deallocate(zkrs)
deallocate(psiperr)
deallocate(zkr)
!
! From allocate_cond_2
!
deallocate(zkrr)
deallocate(newbg)
deallocate(fun0)

2
PWCOND/gep_gep.f90
View File

@ -17,7 +17,7 @@ implicit none
alpha(:), beta(:)
do i=1, n
amt(i,i) = amt(i,i)-delgep
amt(i,i) = amt(i,i)+delgep
bmt(i,i) = bmt(i,i)+delgep
enddo

6
PWCOND/gep_x.f90
View File

@ -11,6 +11,7 @@ subroutine gep_x(n, amt, bmt, eigen, veigen)
! It solves GEP: A X = lambda B X using LAPACK routines
!
USE kinds, only : DP
USE cond, only : delgep
implicit none
integer :: i, n, info, lwork
@ -32,8 +33,8 @@ subroutine gep_x(n, amt, bmt, eigen, veigen)
! pc_ifc, ibmsp and origin. If you have problems, try llapack=.false.
!
do i=1,n
amt(i,i)=amt(i,i)+5.d-10
bmt(i,i)=bmt(i,i)+5.d-10
amt(i,i)=amt(i,i)+delgep
bmt(i,i)=bmt(i,i)+delgep
enddo
call ZGGEV('N', 'V', n, amt, n, bmt, n, alpha, beta, veigen, n, veigen, &
@ -51,6 +52,7 @@ subroutine gep_x(n, amt, bmt, eigen, veigen)
!
do i=1, n
eigen(i)=alpha(i)/beta(i)
! write(6,'(i5, 2f40.20)') i, dreal(eigen(i)), aimag(eigen(i))
enddo
deallocate(work)

562
PWCOND/init_cond.f90
View File

@ -1,306 +1,167 @@
!
! Copyright (C) 2003 A. Smogunov
! 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 init_cond
subroutine init_cond (nregion, flag)
!
! This subroutine sets up some variables of PWCOND
! This subroutine computes and allocates the z mesh and
! the local potential for the left/right leads or for the
! scattering region
!
USE ions_base, ONLY : nat, ityp, ntyp => nsp, atm, tau
USE io_global, ONLY : stdout
USE pwcom
USE noncollin_module, ONLY : noncolin
USE uspp_param, ONLY : dion, nbeta, lll, tvanp
USE cond
implicit none
integer :: nt, ib, ir, is, na, iorb, iorb1, m, ih, ih1, ioins, &
iocros, k, ipol, apol
real(kind=DP) :: ledge, redge, ra, zorb1, zorb2
real(kind=DP), parameter :: eps=1.d-8
! input:
! nregion - number of regions to divide the unit cell
! flag - 'l'/'s'/'r'/'t' if the unit cell containes
! the left lead/scat. reg./right lead/all of them
!
use pwcom
use io_files, only : save_file
USE noncollin_module, ONLY : noncolin, npol
USE uspp_param, ONLY : nbrx, nbeta, lll, betar, tvanp
use atom, only: mesh, r
USE ions_base, ONLY : atm, nat, ityp, ntyp => nsp, tau
use cond
!
! We order the atomic orbitals
!
implicit none
!
! Left tip
!
iocros=0
ioins=nocrosl
do na=1, nat
nt=ityp(na)
ih=0
do ib=1, nbeta(nt)
ledge=tau(3,na)-rsph(ib,nt)
redge=tau(3,na)+rsph(ib,nt)
if (tau(3,na).gt.bdl1+eps.and.tau(3,na).le.bdl2+eps) then
if (ledge.gt.bdl1.and.redge.le.bdl2) then
do m=1,2*lll(ib,nt)+1
ioins=ioins+1
ih=ih+1
natih(ioins,1)=na
natih(ioins,2)=ih
itnew(ioins)=nt
nbnew(ioins)=ib
ls(ioins)=lll(ib,nt)
mnew(ioins)=m
do ipol=1,3
taunew(ipol,ioins)=tau(ipol,na)
enddo
enddo
else
do m=1,2*lll(ib,nt)+1
iocros=iocros+1
ih=ih+1
natih(iocros,1)=na
natih(iocros,2)=ih
itnew(iocros)=nt
nbnew(iocros)=ib
ls(iocros)=lll(ib,nt)
mnew(iocros)=m
natih(iocros+noinsl+nocrosl,1)=na
natih(iocros+noinsl+nocrosl,2)=ih
itnew(iocros+noinsl+nocrosl)=nt
nbnew(iocros+noinsl+nocrosl)=ib
ls(iocros+noinsl+nocrosl)=lll(ib,nt)
mnew(iocros+noinsl+nocrosl)=m
do ipol=1,2
taunew(ipol,iocros)=tau(ipol,na)
taunew(ipol,iocros+noinsl+nocrosl)=tau(ipol,na)
enddo
if (ledge.le.bdl1) then
taunew(3,iocros)=tau(3,na)
taunew(3,iocros+noinsl+nocrosl)=tau(3,na)+bdl2-bdl1
else
taunew(3,iocros)=tau(3,na)-(bdl2-bdl1)
taunew(3,iocros+noinsl+nocrosl)=tau(3,na)
endif
enddo
endif
endif
enddo
enddo
character(len=1) :: flag
integer :: nregion, iregion, nrztot, iz, naux, k, mmax, &
nt, ib, ir, nrz1, info, na
integer :: ipol1, ipol2
real(kind=DP), parameter :: epsbeta=1.d-4, eps=1.d-8
real(kind=DP) :: zlen, dz1, dz2, bd1, bd2, bd3, bd4, bd5, &
bmax, ledge, redge
real(kind=DP), allocatable :: ztot(:), rsph(:,:), dwid(:), &
nrzreg(:)
complex(kind=DP), allocatable :: vppottot(:,:,:,:)
character(len=14) :: extension
if (ikind.ne.0) then
!
! If the scattering problem
!
ioins=2*nocrosl+noinsl
do na=1, nat
nt=ityp(na)
ih=0
do ib=1, nbeta(nt)
ledge=tau(3,na)-rsph(ib,nt)
redge=tau(3,na)+rsph(ib,nt)
if (ledge.gt.bdl2.and.redge.le.bdr1) then
do m=1,2*lll(ib,nt)+1
ioins=ioins+1
ih=ih+1
natih(ioins,1)=na
natih(ioins,2)=ih
itnew(ioins)=nt
nbnew(ioins)=ib
ls(ioins)=lll(ib,nt)
mnew(ioins)=m
do ipol=1,3
taunew(ipol,ioins)=tau(ipol,na)
enddo
enddo
endif
enddo
enddo
iocros=2*nocrosl+noinsl+noinss
ioins=iocros+nocrosr
nrx = nr1s
nry = nr2s
nrztot = nr3s
if(nrztot/2*2.eq.nrztot) nrztot = nrztot+1
zlen = at(3,3)
dz1 = zlen/nrztot
sarea = abs(at(1,1)*at(2,2)-at(2,1)*at(1,2))*alat**2
if (abs(ecut2d).le.eps) ecut2d = ecutwfc
if (ikind.eq.1) then
!
! If tips are identical
!
do ib=1, nocrosr
iocros=iocros+1
natih(iocros,1)=natih(ib,1)
natih(iocros,2)=natih(ib,2)
itnew(iocros)=itnew(ib)
nbnew(iocros)=nbnew(ib)
ls(iocros)=ls(ib)
mnew(iocros)=mnew(ib)
do ipol=1,2
taunew(ipol,iocros)=taunew(ipol,ib)
enddo
taunew(3,iocros)=taunew(3,ib)+(bdr1-bdl1)
enddo
else
!
! Two tips are different
!
do na=1, nat
nt=ityp(na)
ih=0
do ib=1, nbeta(nt)
ledge=tau(3,na)-rsph(ib,nt)
redge=tau(3,na)+rsph(ib,nt)
if (ledge.gt.bdr1.and.redge.le.bdr2) then
do m=1,2*lll(ib,nt)+1
ioins=ioins+1
ih=ih+1
natih(ioins,1)=na
natih(ioins,2)=ih
itnew(ioins)=nt
nbnew(ioins)=ib
ls(ioins)=lll(ib,nt)
mnew(ioins)=m
do ipol=1,3
taunew(ipol,ioins)=tau(ipol,na)
enddo
enddo
endif
if (ledge.le.bdr1.and.redge.gt.bdr1) then
do m=1,2*lll(ib,nt)+1
iocros=iocros+1
ih=ih+1
natih(iocros,1)=na
natih(iocros,2)=ih
itnew(iocros)=nt
nbnew(iocros)=ib
ls(iocros)=lll(ib,nt)
mnew(iocros)=m
natih(iocros+noinsr+nocrosr,1)=na
natih(iocros+noinsr+nocrosr,2)=ih
itnew(iocros+noinsr+nocrosr)=nt
nbnew(iocros+noinsr+nocrosr)=ib
ls(iocros+noinsr+nocrosr)=lll(ib,nt)
mnew(iocros+noinsr+nocrosr)=m
do ipol=1,2
taunew(ipol,iocros)=tau(ipol,na)
taunew(ipol,iocros+noinsr+nocrosr)=tau(ipol,na)
enddo
taunew(3,iocros)=tau(3,na)
taunew(3,iocros+noinsr+nocrosr)=tau(3,na)+ &
(bdr2-bdr1)
enddo
endif
enddo
enddo
allocate ( ztot(nrztot+1) )
allocate ( rsph(nbrx, npsx) )
allocate ( dwid(5) )
allocate ( nrzreg(4) )
bd1 = 0.d0
bd3 = zlen
bd4 = zlen
bd5 = zlen
if(flag.eq.'l') then
bd2 = bdl
elseif(flag.eq.'s') then
bd2 = bds
elseif(flag.eq.'r') then
bd2 = bdr
elseif(flag.eq.'t') then
bd2 = bdl
if(nregion.gt.1) then
bd3 = bds
endif
if(nregion.gt.2) then
bd4 = bdr
endif
endif
if(bd2.le.1.d-6) bd2 = zlen
if(bd3.le.1.d-6) bd3 = zlen
if(bd4.le.1.d-6) bd4 = zlen
!
! To form zpseu
!
zpseu=0.d0
if (noncolin) zpseu_nc=0.d0
dwid(1) = bd1
dwid(2) = bd2
dwid(3) = bd3
dwid(4) = bd4
dwid(5) = bd5
iocros=nocrosr
if (ikind.eq.0) iocros=0
do iorb=nocrosl+1, norb-iocros
ib=nbnew(iorb)
nt=itnew(iorb)
if(tvanp(nt).or.lspinorb) then
na=natih(iorb,1)
ih=natih(iorb,2)
do iorb1=nocrosl+1, norb-iocros
if (na.eq.natih(iorb1,1)) then
ih1=natih(iorb1,2)
do is=1, nspin
if (noncolin) then
zpseu_nc(iorb,iorb1,is)=deeq_nc(ih,ih1,na,is)
else
zpseu(iorb,iorb1,is)=deeq(ih,ih1,na,is)
endif
enddo
endif
enddo
else
do is=1, nspin
if (noncolin) then
zpseu_nc(iorb,iorb,is)=dion(ib,ib,nt)
else
zpseu(iorb,iorb,is)=dion(ib,ib,nt)
endif
enddo
nrz1 = 0
do iz = 2, 5
naux=(dwid(iz)+dz1*0.5d0)/dz1-nrz1
nrzreg(iz-1) = naux
if (naux.gt.0) then
dz2=(dwid(iz)-dwid(iz-1))/naux
do k=1, naux
ztot(nrz1+k)=dwid(iz-1)+(k-1)*dz2
enddo
nrz1 = nrz1 + naux
endif
enddo
enddo
if(nrz1.ne.nrztot) CALL errore ('in init_cond','wrong nrztot',info)
ztot(nrztot+1) = zlen
do iorb=1, nocrosl
do iorb1=1, nocrosl
do is=1, nspin
if (noncolin) then
zpseu_nc(iorb,iorb1,is)= &
zpseu_nc(nocrosl+noinsl+iorb,nocrosl+noinsl+iorb1,is)
else
zpseu(iorb,iorb1,is)= &
zpseu(nocrosl+noinsl+iorb,nocrosl+noinsl+iorb1,is)
endif
enddo
enddo
enddo
allocate (vppottot(nrztot, nrx*nry, npol, npol))
call poten(vppottot,nrztot,ztot)
iocros=0
if (ikind.eq.2) iocros=2*nocrosl+noinsl+noinss
if (ikind.ne.0) then
do iorb=1, nocrosr
do iorb1=1, nocrosr
do is=1, nspin
if (noncolin) then
zpseu_nc(norb-nocrosr+iorb,norb-nocrosr+iorb1,is)= &
zpseu_nc(iocros+iorb,iocros+iorb1,is)
else
zpseu(norb-nocrosr+iorb,norb-nocrosr+iorb1,is)= &
zpseu(iocros+iorb,iocros+iorb1,is)
endif
enddo
enddo
enddo
endif
!
! to form the array containing the information does the orbital
! cross the given slab or not.
! to determine radii of nonlocal spheres
!
do iorb=1, norb
ra=rsph(nbnew(iorb),itnew(iorb))
zorb1=taunew(3,iorb)-ra
zorb2=taunew(3,iorb)+ra
do k=1, nrz
if (zorb1.gt.(z(k+1)-eps).or.zorb2.lt.(z(k)+eps)) then
cross(iorb,k)=0
else
cross(iorb,k)=1
endif
mmax = 0
do nt=1, ntyp
do ib=1, nbeta(nt)
mmax = max(mmax, lll(ib, nt))
bmax=0.d0
do ir=2, mesh(nt)
bmax=max(bmax, betar(ir,ib,nt)/r(ir,nt))
enddo
ir=mesh(nt)
do while (abs(betar(ir,ib,nt)/r(ir,nt)).le.epsbeta*bmax)
ir=ir-1
enddo
rsph(ib,nt)=r(ir,nt)/alat
enddo
enddo
if (mmax.gt.3) call errore ('allocate','for l>3 -orbitals &
& the method is not implemented',1)
!
! Some output by PWCOND
! We set up the radii of US spheres to be the same (to avoid
! the problem with the spheres crossing or not the boundaries)
!
WRITE( stdout,'(9x, ''PWCOND calculation is now starting...'')')
if(ikind.eq.0) then
WRITE( stdout,'(/,2x,''CBS calculation (ikind=0)'')')
WRITE( stdout,'(/,2x,''left boundary is '',f9.4)') bdl1
WRITE( stdout,'(2x,''right boundary is '',f9.4)') bdl2
elseif(ikind.eq.1) then
WRITE( stdout,'(/,2x,''T calculation with identical tips (ikind=1)'')')
WRITE( stdout,'(/,2x,''left boundary of the tip is '',f9.4)') bdl1
WRITE( stdout,'(2x,''right boundary of the tip is '',f9.4)') bdl2
WRITE( stdout,'(2x,''the end of the scatt. region is '',f9.4)') bdr1
elseif(ikind.eq.2) then
WRITE( stdout,'(/,2x,''T calculation with different tips (ikind=2)'')')
WRITE( stdout,'(/,2x,''left boundary of the left tip is '',f9.4)') bdl1
WRITE( stdout,'(2x,''right boundary of the left tip is '',f9.4)') bdl2
WRITE( stdout,'(2x,''left boundary of the right tip is '',f9.4)') bdr1
WRITE( stdout,'(2x,''right boundary of the right tip is '',f9.4)') bdr2
do nt=1, ntyp
if (tvanp(nt)) then
bmax=0.d0
do ib=1, nbeta(nt)
bmax=max(bmax, rsph(ib,nt))
enddo
do ib=1, nbeta(nt)
rsph(ib,nt)=bmax
enddo
endif
enddo
!
! Move all atoms into the unit cell
!
do na=1, nat
if(tau(3,na).gt.zlen) tau(3,na)=tau(3,na)-zlen
if(tau(3,na).le.0) tau(3,na)=tau(3,na)+zlen
enddo
!----------------
! Some output
write(6,*)
if(flag.eq.'l') then
write(6,'(''===== INPUT FILE containing the left lead ====='')')
elseif(flag.eq.'s') then
write(6,'(''===== INPUT FILE containing the scat. region ====='')')
elseif(flag.eq.'r') then
write(6,'(''===== INPUT FILE containing the right lead ====='')')
elseif(flag.eq.'t') then
write(6,'(''===== INPUT FILE containing all the regions ====='')')
endif
WRITE( stdout,'(/,5x,''GEOMETRY:'')')
WRITE( stdout, 100) alat, omega, sarea, zl, nat, ntyp
write(6,'(/,5x,''GEOMETRY:'')')
write (6, 100) alat, omega, sarea, zlen, nat, ntyp
100 format (/,5x, &
& 'lattice parameter (a_0) = ',f12.4,' a.u.',/,5x, &
& 'the volume = ',f12.4,' (a.u.)^3',/,5x, &
@ -309,19 +170,18 @@ subroutine init_cond
& 'number of atoms/cell = ',i12,/,5x, &
& 'number of atomic types = ',i12,/,5x)
WRITE( stdout,'(5x,''crystal axes: (cart. coord. in units of a_0)'',/, &
write(6,'(5x,''crystal axes: (cart. coord. in units of a_0)'',/, &
& 3(15x,''a('',i1,'') = ('',3f8.4,'' ) '',/ ) )') &
& ( apol, (at(ipol,apol), ipol=1,3), apol=1,3)
& ( na, (at(nt,na), nt=1,3), na=1,3)
WRITE( stdout,'(/,3x,''Cartesian axes'')')
WRITE( stdout, '(/,5x,''site n. atom '', &
write(6,'(/,3x,''Cartesian axes'')')
write(6, '(/,5x,''site n. atom '', &
& '' positions (a_0 units)'')')
WRITE( stdout, '(7x,i3,8x,a6,'' tau('',i2,'')=('',3f8.4,'' )'')') &
write(6, '(7x,i4,8x,a6,'' tau('',i3,'')=('',3f8.4,'' )'')') &
& ( na,atm(ityp(na)),na, &
& ( tau(ipol,na),ipol=1,3),na=1,nat )
WRITE( stdout, 300) nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, &
nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrx, nry, nrz, nz1
& ( tau(nt,na),nt=1,3),na=1,nat )
write (6, 300) nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, &
nr1, nr2, nr3, nrx1, nrx2, nrx3
300 format (/,5x, &
& 'nr1s = ',i12,/,5x, &
& 'nr2s = ',i12,/,5x, &
@ -334,63 +194,103 @@ subroutine init_cond
& 'nr3 = ',i12,/,5x, &
& 'nrx1 = ',i12,/,5x, &
& 'nrx2 = ',i12,/,5x, &
& 'nrx3 = ',i12,/,5x, &
& 'nrx = ',i12,/,5x, &
& 'nry = ',i12,/,5x, &
& 'nrz = ',i12,/,5x, &
& 'nz1 = ',i12,/,5x)
& 'nrx3 = ',i12,/,5x)
WRITE( stdout,*) '_______________________________'
WRITE( stdout,*) ' Radii of nonlocal spheres: '
WRITE( stdout, '(/,5x,''type ibeta ang. mom.'', &
write(6,*) '_______________________________'
write(6,*) ' Radii of nonlocal spheres: '
write(6, '(/,5x,''type ibeta ang. mom.'', &
& '' radius (a_0 units)'')')
WRITE( stdout, '(7x,a6,3x,i3,7x,i3,14x,f12.4)') &
write(6, '(7x,a6,3x,i3,7x,i3,14x,f12.4)') &
& ( ( atm(nt), ib, lll(ib,nt), rsph(ib,nt), &
& ib=1,nbeta(nt) ), nt=1,ntyp)
WRITE( stdout, 200) norb, norbf, nocrosl, noinsl, noinss, nocrosr, &
noinsr
200 format (/,5x, &
& 'norb = ',i12,/,5x, &
& 'norbf = ',i12,/,5x, &
& 'nocrosl = ',i12,/,5x, &
& 'noinsl = ',i12,/,5x, &
& 'noinss = ',i12,/,5x, &
& 'nocrosr = ',i12,/,5x, &
& 'noinsr = ',i12,/,5x)
!-----------------------------
WRITE( stdout, '(5x,''iorb type ibeta ang. mom.'',3x, &
& ''m position (a_0)'')')
WRITE( stdout,'(5x,i3,4x,a6,i3,6x,i3,6x,i3,'' taunew('', &
& i3,'')=('',3f8.4,'')'')') &
& ( iorb,atm(itnew(iorb)), nbnew(iorb), ls(iorb), &
& mnew(iorb), iorb, &
& (taunew(ipol,iorb),ipol=1,3), iorb=1, norb )
WRITE( stdout, 301) ef*rytoev, energy0, denergy, nenergy, ecut2d, &
ewind, epsproj
301 format (/,5x, &
& 'Fermi energy = ',1pe15.5,/,5x, &
& 'energy0 = ',1pe15.1,/,5x, &
& 'denergy = ',1pe15.1,/,5x, &
& 'nenergy = ',i10,/,5x, &
& 'ecut2d = ',1pe15.1,/,5x, &
& 'ewind = ',1pe15.1,/,5x, &
& 'epsproj = ',1pe15.1,/,5x)
if(nspin.ne.1) then
WRITE( stdout,'(/,9x, ''Calculations are done for: '')')
WRITE( stdout,'(2x, ''spin index = '', i6)') iofspin
if(flag.eq.'l') then
nrzl = nrzreg(1)
allocate( vppotl(nrzl, nrx * nry, npol, npol) )
allocate( zl(nrzl+1) )
call potz_split(vppottot,ztot,vppotl,zl,nrztot,nrzl,nrx*nry,npol,0)
call init_orbitals(zlen,bd1,bd2,zl,nrzl,rsph,1)
efl = ef
elseif(flag.eq.'s') then
nrzs = nrzreg(1)
allocate( vppots(nrzs, nrx * nry, npol, npol) )
allocate( zs(nrzs+1) )
call potz_split(vppottot,ztot,vppots,zs,nrztot,nrzs,nrx*nry,npol,0)
call init_orbitals(zlen,bd1,bd2,zs,nrzs,rsph,2)
efs = ef
elseif(flag.eq.'r') then
nrzr = nrzreg(1)
allocate( vppotr(nrzr, nrx * nry, npol, npol) )
allocate( zr(nrzr+1) )
call potz_split(vppottot,ztot,vppotr,zr,nrztot,nrzr,nrx*nry,npol,0)
call init_orbitals(zlen,bd1,bd2,zr,nrzr,rsph,3)
efr = ef
elseif(flag.eq.'t') then
nrzl = nrzreg(1)
allocate( vppotl(nrzl, nrx * nry, npol, npol) )
allocate( zl(nrzl+1) )
call potz_split(vppottot,ztot,vppotl,zl,nrztot,nrzl,nrx*nry,npol,0)
call init_orbitals(zlen,bd1,bd2,zl,nrzl,rsph,1)
if(nregion.gt.1) then
nrzs = nrzreg(2)
allocate( vppots(nrzs, nrx * nry, npol, npol) )
allocate( zs(nrzs+1) )
call potz_split(vppottot,ztot,vppots,zs,nrztot,nrzs,nrx*nry, &
npol,nrzl)
call init_orbitals(zlen,bd2,bd3,zs,nrzs,rsph,2)
endif
if(nregion.gt.2) then
nrzr = nrzreg(3)
allocate( vppotr(nrzr, nrx * nry, npol, npol) )
allocate( zr(nrzr+1) )
call potz_split(vppottot,ztot,vppotr,zr,nrztot,nrzr,nrx*nry, &
npol,nrzl+nrzs)
call init_orbitals(zlen,bd3,bd4,zr,nrzr,rsph,3)
endif
efl = ef
efs = ef
efr = ef
endif
if(norb.le.80) then
WRITE( stdout,'(4x,''k slab'',3x,'' z(k) z(k+1)'', &
& 5x,''crossing(iorb=1,norb)'')')
do k=1, nrz
WRITE( stdout,'(2x,i3,2x,3f7.4,3x,80i1)') &
k,z(k),z(k+1),z(k+1)-z(k),(cross(iorb,k),iorb=1,norb)
enddo
endif
deallocate(dwid)
deallocate (ztot)
deallocate (rsph)
deallocate (nrzreg)
deallocate (vppottot)
return
end subroutine init_cond
subroutine potz_split(vppottot,ztot,vppot,z,nrztot,nrz,nrxy,npol,iz0)
!
! vppottot and ztot --> vppot and z
!
use kinds, only : DP
implicit none
integer :: nrztot, nrz, nrxy, npol, iz0, iz, ixy, ipol1, ipol2
real(kind=DP) :: ztot(nrztot+1), z(nrz+1), zinit
complex(kind=DP) :: vppottot (nrztot, nrxy, npol, npol), &
vppot (nrz, nrxy, npol, npol)
do iz = 1, nrz
do ixy = 1, nrxy
do ipol1 = 1, npol
do ipol2 = 1, npol
vppot(iz,ixy,ipol1,ipol2) = vppottot(iz0+iz,ixy,ipol1,ipol2)
enddo
enddo
enddo
z(1+iz) = ztot(iz0+1+iz) - zinit
enddo
zinit = ztot(iz0+1)
do iz = 1, nrz+1
z(iz) = ztot(iz0+iz) - zinit
enddo
return
end subroutine init_cond
end subroutine potz_split

28
PWCOND/init_gper.f90
View File

@ -21,13 +21,13 @@ subroutine init_gper(ik)
real(kind=DP), parameter :: eps=1.d-8
real(kind=DP), allocatable :: gnorm2(:)
allocate( gnorm2( nr1s * nr2s ) )
allocate( nshell( nr1s, nr2s ) )
allocate( gnorm2( nrx * nry ) )
allocate( nshell( nrx, nry ) )
!
! Compute ngper, ngpsh, and (i,j) --> shell
!
do i=1, nr1s
do j=1, nr2s
do i=1, nrx
do j=1, nry
nshell(i,j)=0
enddo
enddo
@ -36,12 +36,12 @@ subroutine init_gper(ik)
ngpsh=1
gnorm2(1)=(xyk(1,ik)**2+xyk(2,ik)**2)*tpiba2
nshell(1,1)=1
do i=1, nr1s
do i=1, nrx
il=i-1
if (il.gt.nr1s/2) il=il-nr1s
do j=1, nr2s
if (il.gt.nrx/2) il=il-nrx
do j=1, nry
jl=j-1
if (jl.gt.nr2s/2) jl=jl-nr2s
if (jl.gt.nry/2) jl=jl-nry
norm2=(((il+xyk(1,ik))*bg(1,1)+(jl+xyk(2,ik))*bg(1,2))**2+ &
((il+xyk(1,ik))*bg(2,1)+(jl+xyk(2,ik))*bg(2,2))**2)* &
tpiba2
@ -80,8 +80,8 @@ subroutine init_gper(ik)
do i=1, ngpsh
ninsh(i)=0
enddo
do i=1, nr1s
do j=1, nr2s
do i=1, nrx
do j=1, nry
if (nshell(i,j).ne.0) then
do k=nshell(i,j)+1, ngpsh
ninsh(k)=ninsh(k)+1
@ -92,12 +92,12 @@ subroutine init_gper(ik)
!
! To form g
!
do i=1, nr1s
do i=1, nrx
il=i-1
if (il.gt.nr1s/2) il=il-nr1s
do j=1, nr2s
if (il.gt.nrx/2) il=il-nrx
do j=1, nry
jl=j-1
if (jl.gt.nr2s/2) jl=jl-nr2s
if (jl.gt.nry/2) jl=jl-nry
if (nshell(i,j).ne.0) then
ninsh(nshell(i,j))=ninsh(nshell(i,j))+1
igper=ninsh(nshell(i,j))

362
PWCOND/init_orbitals.f90 Normal file
View File

@ -0,0 +1,362 @@
!
! 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 init_orbitals (zlen, bd1, bd2, z, nrz, rsph, lsr)
!
! Calculates and allocates some variables describing the nonlocal
! orbitals
!
! input:
! zlen - the length of the unit cell in the z direction
! bd1, bd2 - two boundaries of the region under interest
! z(nrz) - mesh in the z direction
! rsph - radii of the spheres
! lsr - 1/2/3 if the region is left lead/scat. reg./right lead
!
use cond
use lsda_mod, only: nspin
use noncollin_module, only : noncolin
use spin_orb, only: lspinorb
use ions_base, only : atm, nat, ityp, tau
use uspp_param, only : nbrx, nbeta, lll, betar, tvanp, dion
use uspp, only : deeq, deeq_nc, qq, qq_so
use atom, only : r, rab
implicit none
integer :: noins, lnocros, rnocros, nocros, norb, na, nt, ih, ih1,&
ioins, ilocros, irocros, orbin, orbfin, ib, lsr, nrz, &
m, k, ipol, iorb, iorb1, is
integer, allocatable :: orbind(:,:), tblm(:,:), cros(:,:), natih(:,:)
real(kind=DP), parameter :: eps=1.d-8
real(kind=DP) :: ledge, redge, ledgel, redgel, ledger, redger, &
bd1, bd2, zlen, z(nrz+1), rsph(nbrx, npsx)
real(kind=DP), allocatable :: taunew(:,:), zpseu(:,:,:)
complex(kind=DP), allocatable :: zpseu_nc(:,:,:,:)
allocate ( orbind(nat,nbrx) )
orbind = -1
!---------------------
! Calculate number of crossing and inside-lying orbitals
!
noins = 0
lnocros = 0
rnocros = 0
do na = 1, nat
nt = ityp(na)
do ib = 1, nbeta(nt)
ledge = tau(3,na)-rsph(ib,nt)
ledgel = ledge-zlen
ledger = ledge+zlen
redge = tau(3,na)+rsph(ib,nt)
redgel = redge-zlen
redger = redge+zlen
if (ledge.le.bd1.and.redge.gt.bd2) &
call errore ('init_cond','Too big atomic spheres',1)
if (ledge.gt.bd1.and.redge.le.bd2) then
noins = noins+2*lll(ib,nt)+1
orbind(na,ib) = 0
elseif(ledge.le.bd1.and.redge.gt.bd1) then
lnocros = lnocros+2*lll(ib,nt)+1
orbind(na,ib) = 1
if(ledger.le.bd2.and.redger.gt.bd2) then
rnocros = rnocros+2*lll(ib,nt)+1
orbind(na,ib) = 2
endif
elseif(ledger.le.bd2.and.redger.gt.bd2) then
rnocros = rnocros+2*lll(ib,nt)+1
orbind(na,ib) = 3
elseif(ledge.le.bd2.and.redge.gt.bd2) then
rnocros = rnocros+2*lll(ib,nt)+1
orbind(na,ib) = 4
if(ledgel.le.bd1.and.redgel.gt.bd1) then
lnocros = lnocros+2*lll(ib,nt)+1
orbind(na,ib) = 5
endif
elseif(ledgel.le.bd1.and.redgel.gt.bd1) then
lnocros = lnocros+2*lll(ib,nt)+1
orbind(na,ib) = 6
endif
enddo
enddo
norb = noins + lnocros + rnocros
nocros = (lnocros + rnocros)/2
!------------------------------------
!-----------------------------
! Formation of some orbital arrays
!
allocate( taunew(4,norb) )
allocate( tblm(4,norb) )
allocate( natih(2,norb) )
allocate( cros(norb, nrz) )
if (noncolin) then
allocate(zpseu_nc(2, norb, norb, nspin))
else
allocate( zpseu(2, norb, norb) )
endif
ilocros = 0
ioins = ilocros + lnocros
irocros = ioins + noins
do na = 1, nat
nt = ityp(na)
ih = 0
do ib = 1, nbeta(nt)
do m = 1,2*lll(ib,nt) + 1
ih = ih+1
if(orbind(na,ib).eq.0) then
ioins = ioins+1
natih(1,ioins)=na
natih(2,ioins)=ih
tblm(1,ioins) = nt
tblm(2,ioins) = ib
tblm(3,ioins) = lll(ib,nt)
tblm(4,ioins) = m
do ipol = 1, 3
taunew(ipol,ioins)=tau(ipol,na)
enddo
taunew(4,ioins) = rsph(ib,nt)
endif
if(orbind(na,ib).eq.1.or.orbind(na,ib).eq.2) then
ilocros = ilocros + 1
natih(1,ilocros)=na
natih(2,ilocros)=ih
tblm(1,ilocros) = nt
tblm(2,ilocros) = ib
tblm(3,ilocros) = lll(ib,nt)
tblm(4,ilocros) = m
do ipol = 1, 3
taunew(ipol,ilocros)=tau(ipol,na)
enddo
taunew(4,ilocros) = rsph(ib,nt)
endif
if(orbind(na,ib).eq.2.or.orbind(na,ib).eq.3) then
irocros = irocros + 1
natih(1,irocros)=na
natih(2,irocros)=ih
tblm(1,irocros) = nt
tblm(2,irocros) = ib
tblm(3,irocros) = lll(ib,nt)
tblm(4,irocros) = m
do ipol = 1, 2
taunew(ipol,irocros)=tau(ipol,na)
enddo
taunew(3,irocros) = tau(3,na) + zlen
taunew(4,irocros) = rsph(ib,nt)
endif
if(orbind(na,ib).eq.4.or.orbind(na,ib).eq.5) then
irocros = irocros + 1
natih(1,irocros)=na
natih(2,irocros)=ih
tblm(1,irocros) = nt
tblm(2,irocros) = ib
tblm(3,irocros) = lll(ib,nt)
tblm(4,irocros) = m
do ipol = 1, 3
taunew(ipol,irocros)=tau(ipol,na)
enddo
taunew(4,irocros) = rsph(ib,nt)
endif
if(orbind(na,ib).eq.5.or.orbind(na,ib).eq.6) then
ilocros = ilocros + 1
natih(1,ilocros)=na
natih(2,ilocros)=ih
tblm(1,ilocros) = nt
tblm(2,ilocros) = ib
tblm(3,ilocros) = lll(ib,nt)
tblm(4,ilocros) = m
do ipol = 1, 2
taunew(ipol,ilocros)=tau(ipol,na)
enddo
taunew(3,ilocros) = tau(3,na) - zlen
taunew(4,ilocros) = rsph(ib,nt)
endif
enddo
enddo
enddo
do iorb = 1, norb
taunew(3,iorb) = taunew(3,iorb) - bd1
enddo
!--------------------------
!-------------------------
! to form the array containing the information does the orbital
! cross the given slab or not.
!
do iorb=1, norb
ledge = taunew(3,iorb)-taunew(4,iorb)
redge = taunew(3,iorb)+taunew(4,iorb)
do k=1, nrz
if (ledge.gt.z(k+1).or.redge.lt.z(k)) then
cros(iorb,k)=0
else
cros(iorb,k)=1
endif
enddo
enddo
!----------------------------
!----------------------------
! To form zpseu
!
if (noncolin) then
zpseu_nc=0.d0
else
zpseu = 0.d0
endif
orbin = 1
orbfin = lnocros+noins
do k = 1, 2
do iorb = orbin, orbfin
nt = tblm(1,iorb)
ib = tblm(2,iorb)
if(tvanp(nt).or.lspinorb) then
na = natih(1,iorb)
ih = natih(2,iorb)
do iorb1 = orbin, orbfin
if (na.eq.natih(1,iorb1)) then
ih1 = natih(2,iorb1)
if (noncolin) then
do is=1, nspin
if(lspinorb) then
zpseu_nc(1,iorb,iorb1,is)=deeq_nc(ih,ih1,na,is)
zpseu_nc(2,iorb,iorb1,is)=qq_so(ih,ih1,is,nt)
else
zpseu_nc(1,iorb,iorb1,is)=deeq_nc(ih,ih1,na,is)
zpseu_nc(2,iorb,iorb1,is)=qq(ih,ih1,nt)
endif
enddo
else
zpseu(1,iorb,iorb1)=deeq(ih,ih1,na,iofspin)
zpseu(2,iorb,iorb1) = qq(ih,ih1,nt)
endif
endif
enddo
else
if (noncolin) then
do is = 1, nspin
zpseu_nc(1,iorb,iorb,is)=dion(ib,ib,nt)
enddo
else
zpseu(1,iorb,iorb)=dion(ib,ib,nt)
endif
endif
enddo
orbin = lnocros+noins+1
orbfin = norb
enddo
!--------------------------
!--------------------------
! Allocation
!
if(lsr.eq.1) then
norbl = norb
nocrosl = nocros
noinsl = noins
if(ikind.eq.1) then
norbr = norb
nocrosr = nocros
noinsr = noins
endif
allocate( taunewl(4,norbl) )
allocate( tblml(4,norbl) )
allocate( crosl(norbl, nrzl) )
if (noncolin) then
allocate(zpseul_nc(2, norbl, norbl, nspin))
else
allocate( zpseul(2, norbl, norbl) )
endif
taunewl = taunew
tblml = tblm
crosl = cros
if (noncolin) then
zpseul_nc = zpseu_nc
else
zpseul = zpseu
endif
rl = r
rabl = rab
betarl = betar
norbf = norbl
elseif(lsr.eq.2) then
norbs = norb
noinss = noins
allocate( taunews(4,norbs) )
allocate( tblms(4,norbs) )
allocate( cross(norbs, nrzs) )
if (noncolin) then
allocate(zpseus_nc(2, norbs, norbs, nspin))
else
allocate( zpseus(2, norbs, norbs) )
endif
taunews = taunew
tblms = tblm
cross = cros
if (noncolin) then
zpseus_nc = zpseu_nc
else
zpseus = zpseu
endif
rs = r
rabs = rab
betars = betar
norbf = max(norbf,norbs)
elseif(lsr.eq.3) then
norbr = norb
nocrosr = nocros
noinsr = noins
allocate( taunewr(4,norbr) )
allocate( tblmr(4,norbr) )
allocate( crosr(norbr, nrzr) )
if (noncolin) then
allocate(zpseur_nc(2, norbr, norbr, nspin))
else
allocate( zpseur(2, norbr, norbr) )
endif
taunewr = taunew
tblmr = tblm
crosr = cros
if (noncolin) then
zpseur_nc = zpseu_nc
else
zpseur = zpseu
endif
rr = r
rabr = rab
betarr = betar
norbf = max(norbf,norbr)
endif
!---------------------------
deallocate (orbind)
deallocate (taunew)
deallocate (tblm)
deallocate (natih)
deallocate (cros)
if (noncolin) then
deallocate (zpseu_nc)
else
deallocate (zpseu)
endif
return
end subroutine init_orbitals

40
PWCOND/jbloch.f90
View File

@ -9,8 +9,8 @@
! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC).
!
!
subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
vec, kval, c1, c2, sarea, nchan, npol)
subroutine jbloch (nst, n2d, norbf, norb, nocros, kfun, kfund, &
vec, kval, c1, c2, nchan, npol)
!
! This routine computes the current I carrying by the Bloch state
! so that <psi_k|I|psi_k'> = \delta_{kk'} I_k and does some
@ -19,6 +19,7 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
#include "f_defs.h"
USE kinds, only : DP
USE noncollin_module, only : noncolin
USE cond, only : sarea
implicit none
real(kind=DP), parameter :: eps=1.d-4
@ -27,19 +28,18 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
n2d, & ! 2D-dimension
noins, & ! number of interior orbitals
nocros, & ! number of the orbitals crossing the boundary
norbnow, & ! total number of orbitals =noins+2*nocros
norb, & ! total number of orbitals =noins+2*nocros
norbf, & ! max number of orbitals
npol, & ! number of wave-functions components
nst, & ! number of Bloch states =2*(n2d+nocros)
nchan, & ! number of propagating channels
info, i, j, k, n, iorb, n1, il, ir, in
real(kind=DP) :: &
sarea, & ! cross-section
k1, k2
complex(kind=DP) :: &
kfun(n2d, nst), & ! phi_k(d)
kfund(n2d, nst), & ! phi'_k(d)
vec(2*n2d+npol*norbnow, nst), & ! exp. coeff. for phi_k
vec(2*n2d+npol*norb, nst), & ! exp. coeff. for phi_k
kval(nst), & ! k of phi_k
c1(norbf*npol,2*n2d), c2(norbf*npol,norbf*npol), & ! nonlocal integrals
z, z1, ZDOTC
@ -49,7 +49,7 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
complex(kind=DP), allocatable :: kval1(:), vec1(:,:), kcoef(:,:), &
kcuroff(:,:), valj(:,:)
noins = norbnow-2*nocros
noins = norb-2*nocros
!
! To compute the number of channels and ncond(+-,>,<) --> (...)
!
@ -69,7 +69,7 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
allocate( valj( 2*nchan, 2*nchan ) )
allocate( kval1( nst ) )
allocate( kcur( nst ) )
allocate( vec1( (2*n2d+npol*norbnow), nst ) )
allocate( vec1( (2*n2d+npol*norb), nst ) )
kcur=0.d0
kcoef=(0.d0,0.d0)
@ -81,7 +81,7 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
kcoef(iorb,k)=kcoef(iorb,k)+ &
vec(j,ir)*c1(npol*(nocros+noins)+iorb,j)
enddo
do j=1, norbnow*npol
do j=1, norb*npol
kcoef(iorb,k)=kcoef(iorb,k)+ &
vec(2*n2d+j,ir)*c2(npol*(nocros+noins)+iorb,j)
enddo
@ -125,7 +125,7 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
do in=1, 2*nchan
if (ej(in).gt.0.d0) then
ir=ir+1
do n=1, 2*n2d+norbnow*npol
do n=1, 2*n2d+norb*npol
vec1(n,ir)=0.d0
do j=1, 2*nchan
vec1(n,ir)=vec1(n,ir)+vec(n,ncond(j))*valj(j,in)
@ -138,7 +138,7 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
enddo
else
il=il+1
do n=1, 2*n2d+npol*norbnow
do n=1, 2*n2d+npol*norb
vec1(n,il)=0.d0
do j=1, 2*nchan
vec1(n,il)=vec1(n,il)+vec(n,ncond(j))*valj(j,in)
@ -152,27 +152,17 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
endif
enddo
! do j=1,2*nchan
! write(6,'("------------------------------",i5,2f15.6)') j, kval1(j)
! do k=1,n2d
! if (k.eq.n2d/2+1.and.noncolin) write(6,'("-------")')
! write(6,'(i5,f15.7)') k, abs(kfun(k,ncond(j)))
! enddo
! enddo
!-- decaying states
do in=1, nst
if (DIMAG(kval(in)).gt.eps) then
ir=ir+1
kval1(ir)=kval(in)
call DCOPY(2*(2*n2d+npol*norbnow),vec(1,in),1,vec1(1,ir),1)
call DCOPY(2*(2*n2d+npol*norb),vec(1,in),1,vec1(1,ir),1)
endif
if (-DIMAG(kval(in)).gt.eps) then
il=il+1
kval1(il)=kval(in)
call DCOPY(2*(2*n2d+npol*norbnow),vec(1,in),1,vec1(1,il),1)
call DCOPY(2*(2*n2d+npol*norb),vec(1,in),1,vec1(1,il),1)
endif
enddo
@ -181,11 +171,11 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
!
do k=1, nchan
k1=1.d0/abs(kcur(k))
call DSCAL(2*(2*n2d+npol*norbnow),sqrt(k1),vec1(1,k),1)
call DSCAL(2*(2*n2d+npol*norb),sqrt(k1),vec1(1,k),1)
enddo
do k=nst/2+1, nst/2+nchan
k1=1.d0/abs(kcur(k))
call DSCAL(2*(2*n2d+npol*norbnow),sqrt(k1),vec1(1,k),1)
call DSCAL(2*(2*n2d+npol*norb),sqrt(k1),vec1(1,k),1)
enddo
!
@ -221,7 +211,7 @@ subroutine jbloch (nst, n2d, norbf, norbnow, nocros, kfun, kfund, &
kval(k)=kval1(k)
enddo
do k=1, nst
call DCOPY(2*(2*n2d+npol*norbnow),vec1(1,ncond(k)),1,vec(1,k),1)
call DCOPY(2*(2*n2d+npol*norb),vec1(1,ncond(k)),1,vec(1,k),1)
enddo
deallocate(kcoef)

2
PWCOND/kbloch.f90
View File

@ -19,7 +19,7 @@ subroutine kbloch(ntot, val)
integer :: &
ntot, & ! number of Bloch states
in
real(kind=DP) :: d, tpiba, rho, g1, g2, k1, k2
real(kind=DP) :: rho, g1, g2, k1, k2
complex(kind=DP) :: &
val(ntot) ! complex k values

320
PWCOND/local.f90
View File

@ -13,67 +13,183 @@ SUBROUTINE local
!
! This subroutine computes 2D eigenfunctions and eigenvalues for
! the local potential in each slab and performs 2D reduction of
! the plane wave basis set.
! the plane wave basis set (local_1). Using this reduced basis
! set it solves again 2D EV problem (local_2).
!
USE io_global, ONLY : stdout, ionode, ionode_id
USE pwcom
USE constants, ONLY : rytoev
USE io_global, ONLY : stdout
USE noncollin_module, ONLY : npol
USE io_files
USE cond
USE mp, ONLY : mp_barrier, mp_bcast
USE mp_global, ONLY : nproc, me_pool, root_pool
USE parallel_include
IMPLICIT NONE
INTEGER :: i, il, j, jl, ixy, ig, jg, ipol, igper, k, &
ios, index, number, nprob, nteam, nteamnow, &
status, info, kin, kfin, is, js
INTEGER, ALLOCATABLE :: fftxy(:,:)
REAL(kind=DP), PARAMETER :: eps=1.d-6
REAL(kind=DP), ALLOCATABLE :: el(:), gp(:)
COMPLEX(kind=DP), ALLOCATABLE :: amat(:,:), amat1(:,:), ymat(:,:), &
psibase(:,:), psiprob(:,:)
COMPLEX(kind=DP),PARAMETER :: one=(1.d0,0.d0), zero=(0.d0,0.d0)
COMPLEX(kind=DP) :: aij, xfact, ZDOTC
INTEGER :: ig, il, k, kin, kfin
REAL(kind=DP) :: edummy
COMPLEX(kind=DP), ALLOCATABLE :: psibase(:,:)
LOGICAL :: exst
!
! To divide the slabs between CPU
!
CALL start_clock('local')
CALL slabcpu(nrz, nrzp, nkofz, bdl1, bdl2, bdr1, bdr2, z)
!
! If all the information is already contained in the file it reads it.
!
IF (lread_loc) THEN
IF (lread_loc) THEN
CALL seqopn(4,fil_loc,'unformatted',exst)
IF(.NOT.exst) CALL errore ('local','fil_loc not found',1)
READ(4) n2d
READ(4) nrzpl, nrzps, nrzpr
! Allocate variables depending on n2d
CALL allocate_cond_2
READ(4) ((newbg(ig,il), ig=1, ngper*npol), il=1, n2d)
! WRITE( stdout,*) 'ngper, n2d = ', ngper, n2d
READ(4) (((psiper(ig,il,k),ig=1,n2d),il=1,n2d), &
k=1,nrzp)
READ(4) ((zkr(il,k),il=1,n2d),k=1,nrzp)
CALL allocate_cond
READ(4) ((newbg(ig,il), ig=1, ngper*npol), il=1, n2d)
READ(4) (((psiperl(ig,il,k),ig=1,n2d),il=1,n2d), &
k=1,nrzpl)
READ(4) ((zkrl(il,k),il=1,n2d),k=1,nrzpl)
if(ikind.gt.0) then
READ(4) (((psipers(ig,il,k),ig=1,n2d),il=1,n2d), &
k=1,nrzps)
READ(4) ((zkrs(il,k),il=1,n2d),k=1,nrzps)
endif
if(ikind.gt.1) then
READ(4) (((psiperr(ig,il,k),ig=1,n2d),il=1,n2d), &
k=1,nrzpr)
READ(4) ((zkrr(il,k),il=1,n2d),k=1,nrzpr)
endif
CLOSE(unit=4)
RETURN
ENDIF
allocate( psibase( ngper*npol, ngper*npol ) )
psibase = 0.d0
if(ewind.le.100.d0) then
n2d = 0
edummy = earr(1)/rytoev + efl
call local_1(edummy,nrzl,vppotl,n2d,psibase)
if(ikind.gt.0) then
edummy = earr(1)/rytoev + efs
call local_1(edummy,nrzs,vppots,n2d,psibase)
endif
if(ikind.eq.2) then
edummy = earr(1)/rytoev + efr
call local_1(edummy,nrzr,vppotr,n2d,psibase)
endif
else
n2d = ngper*npol
endif
!
! Allocate variables depending on n2d
!
nrzps = 0
nrzpr = 0
call divide(nrzl, kin, kfin)
nrzpl = kfin - kin + 1
if(ikind.gt.0) then
call divide(nrzs, kin, kfin)
nrzps = kfin - kin + 1
endif
if(ikind.gt.1) then
call divide(nrzr, kin, kfin)
nrzpr = kfin - kin + 1
endif
CALL allocate_cond
IF (npol.EQ.2) THEN
WRITE( stdout,*) 'ngper, ngper*npol, n2d = ', ngper, ngper*npol, n2d
ELSE
WRITE( stdout,*) 'ngper, n2d = ', ngper, n2d
ENDIF
!
! Construct components of basis vector set on G_per
!
if(ewind.le.100.d0) then
CALL DCOPY(2*ngper*npol*n2d,psibase,1,newbg,1)
else
newbg = 0.d0
do ig=1, n2d
newbg(ig,ig) = 1.d0
enddo
endif
deallocate( psibase )
call local_2(nrzl,nrzpl,vppotl,psiperl,zkrl)
if(ikind.gt.0) call local_2(nrzs,nrzps,vppots,psipers,zkrs)
if(ikind.gt.1) call local_2(nrzr,nrzpr,vppotr,psiperr,zkrr)
!
! saving the 2D data on the file if lwrite_loc=.t.
!
IF (lwrite_loc) THEN
IF(fil_loc.EQ.' ') CALL errore ('local','fil_loc no name',1)
CALL seqopn(4,fil_loc,'unformatted',exst)
WRITE(4) n2d
WRITE(4) nrzpl, nrzps, nrzpr
WRITE(4) ((newbg(ig,il), ig=1, ngper*npol), il=1, n2d)
WRITE(4) (((psiperl(ig,il,k),ig=1,n2d),il=1,n2d), &
k=1,nrzpl)
WRITE(4) ((zkrl(il,k),il=1,n2d),k=1,nrzpl)
if(ikind.gt.0) then
WRITE(4) (((psipers(ig,il,k),ig=1,n2d),il=1,n2d), &
k=1,nrzps)
WRITE(4) ((zkrs(il,k),il=1,n2d),k=1,nrzps)
endif
if(ikind.gt.1) then
WRITE(4) (((psiperr(ig,il,k),ig=1,n2d),il=1,n2d), &
k=1,nrzpr)
WRITE(4) ((zkrr(il,k),il=1,n2d),k=1,nrzpr)
endif
CLOSE(unit=4)
ENDIF
CALL stop_clock('local')
RETURN
END SUBROUTINE local
!-----------------------------------
subroutine local_1 (edummy, nrz, vppot, n2d, psibase)
USE kinds, only : DP
USE cell_base, ONLY : at, tpiba2
USE noncollin_module, ONLY : npol
USE mp_global, ONLY : nproc, me_pool, root_pool
USE mp, ONLY : mp_barrier, mp_bcast
USE io_global, ONLY : ionode, ionode_id
USE parallel_include
use cond, only : nrx, nry, ngper, gper, ewind, epsproj
IMPLICIT NONE
INTEGER :: nrz, n2d
INTEGER :: i, il, j, jl, ixy, ig, jg, ipol, igper, k, &
ios, index, number, nprob, nteam, nteamnow, &
status, info, kin, kfin, is, js
INTEGER, ALLOCATABLE :: fftxy(:,:)
REAL(kind=DP) :: edummy
REAL(kind=DP), ALLOCATABLE :: el(:), gp(:)
complex(kind=DP) :: psibase(ngper*npol,ngper*npol), &
vppot(nrz,nrx*nry,npol,npol), aij, xfact, ZDOTC
COMPLEX(kind=DP), ALLOCATABLE :: amat(:,:), psiprob(:,:)
COMPLEX(kind=DP),PARAMETER :: one=(1.d0,0.d0), zero=(0.d0,0.d0)
ALLOCATE( gp( 2 ) )
ALLOCATE( el( ngper * npol ) )
ALLOCATE( amat( ngper * npol, ngper * npol ) )
ALLOCATE( psibase( ngper * npol, ngper * npol ) )
ALLOCATE( psiprob( ngper * npol, ngper * npol ) )
ALLOCATE( fftxy(-(nrx-1)/2:nrx/2,-(nry-1)/2:nry/2) )
ALLOCATE( fftxy(-nrx:nrx,-nry:nry) )
!
! To form fftxy correspondence
!
fftxy = 0
DO i=1, nrx
il=i-1
IF (il.GT.nrx/2) il=il-nrx
@ -84,30 +200,19 @@ SUBROUTINE local
ENDDO
ENDDO
!
! to find kin and kfin
!
DO k=1, nrz
IF (z(k).LE.bdl1+eps) kin=k
IF (z(k).LE.bdr2-eps) kfin=k
ENDDO
!
! Starting k and number of CPU
!
nteam=1
kin = 1
kfin = nrz
kin = kin + me_pool
nteam = nproc
!
! set and solve the eigenvalue equation for each slab
!
n2d=0
nprob=0
psibase=(0.d0,0.d0)
DO WHILE(kin.LE.kfin)
amat=(0.d0,0.d0)
DO ig=1, ngper
DO jg=ig, ngper
@ -120,7 +225,7 @@ SUBROUTINE local
DO js=1,npol
amat(ig+(is-1)*ngper,jg+(js-1)*ngper)=vppot(kin,index,is,js)
amat(jg+(js-1)*ngper,ig+(is-1)*ngper)= &
DCONJG(amat(ig+(is-1)*ngper,jg+(js-1)*ngper))
CONJG(amat(ig+(is-1)*ngper,jg+(js-1)*ngper))
ENDDO
ENDDO
ENDIF
@ -132,19 +237,10 @@ SUBROUTINE local
ENDDO
ENDDO
CALL hev_ab(ngper*npol, amat, ngper*npol, el, psiprob, &
-1.d1, eryd+ewind, nprob)
! do is=1,ngper*npol
! write(6,'("------------------------------",i5,f15.7)') is, el(is)
! do ig=1,ngper*npol
! if (ig.eq.ngper+1) write(6,'("----")')
! write(6,'(i5,2f15.7)') ig, psiprob(ig,is)
! enddo
! enddo
! stop
-1.d1, edummy+ewind, nprob)
#ifdef __PARA
IF ( me_pool .NE. root_pool ) THEN
IF ( me_pool.ne.root_pool ) THEN
CALL mpi_send(nprob,1,MPI_INTEGER,0,17, &
MPI_COMM_WORLD,info )
CALL errore ('n2d reduction','info<>0 in send',info)
@ -173,37 +269,70 @@ SUBROUTINE local
CALL gramsh(ngper*npol,nprob,1,nprob,psibase,psiprob,n2d,epsproj)
#endif
kin=kin+nteam
ENDDO
#ifdef __PARA
CALL mp_barrier()
CALL mp_bcast(n2d,ionode_id)
CALL mp_bcast(psibase,ionode_id)
CALL mp_bcast(psibase,ionode_id)
#endif
!
! Allocate variables depending on n2d
!
CALL allocate_cond_2
IF (npol.EQ.2) THEN
WRITE( stdout,*) 'ngper, ngper*npol, n2d = ', ngper, ngper*npol, n2d
ELSE
WRITE( stdout,*) 'ngper, n2d = ', ngper, n2d
ENDIF
!
! Construct components of basis vector set on G_per
!
CALL DCOPY(2*ngper*npol*n2d,psibase,1,newbg,1)
deallocate( gp )
deallocate( el )
deallocate( amat )
deallocate( psiprob )
deallocate( fftxy )
!
! set and solve the eigenvalue equation for each slab
!
return
end subroutine local_1
subroutine local_2(nrz, nrzp, vppot, psiper, zkr)
USE kinds, only : DP
USE cell_base, ONLY : at, tpiba2
USE mp, ONLY : mp_barrier
USE noncollin_module, ONLY : npol
use cond, only : nrx, nry, ngper, n2d, gper, newbg
IMPLICIT NONE
INTEGER :: nrz, nrzp
INTEGER :: i, il, j, jl, ixy, ig, jg, ipol, igper, k, kp, &
ios, info, index, number, kin, kfin, is, js
INTEGER, ALLOCATABLE :: fftxy(:,:)
REAL(kind=DP) :: zkr(n2d,nrzp)
REAL(kind=DP), ALLOCATABLE :: gp(:)
complex(kind=DP) :: psiper(n2d,n2d,nrzp), &
vppot(nrz,nrx*nry,npol,npol), aij, ZDOTC
COMPLEX(kind=DP), ALLOCATABLE :: amat(:,:), amat1(:,:), ymat(:,:)
COMPLEX(kind=DP),PARAMETER :: one=(1.d0,0.d0), zero=(0.d0,0.d0)
allocate( gp( 2 ) )
allocate( fftxy(-nrx:nrx, -nry:nry) )
ALLOCATE( amat( ngper * npol, ngper * npol ) )
ALLOCATE( amat1( n2d, n2d ) )
ALLOCATE( ymat( ngper*npol, n2d ) )
!
! To form fftxy correspondence
!
fftxy(:,:) = 0
do i = 1, nrx
il = i-1
if (il.gt.nrx/2) il=il-nrx
do j = 1, nry
jl = j-1
if (jl.gt.nry/2) jl = jl-nry
fftxy(il,jl) = i+(j-1)*nrx
enddo
enddo
call divide(nrz, kin, kfin)
! for reduced basis set H'_{ab}=e*^i_aH_{ij}e^j_b
DO k=1, nrz
IF(nkofz(k).NE.0) THEN
do k = kin, kfin
kp = k - kin + 1
ymat=(0.d0,0.d0)
!
! First compute y_{ib}=H_{ij}e_{jb}
@ -244,40 +373,24 @@ SUBROUTINE local
! Solving the eigenvalue problem and construction zk
!
info=-1
CALL hev_ab(n2d, amat1, n2d, zkr(1,nkofz(k)), &
psiper(1,1,nkofz(k)), 0.d0, 0.d0, info)
CALL hev_ab(n2d, amat1, n2d, zkr(1,kp), &
psiper(1,1,kp), 0.d0, 0.d0, info)
ENDIF
ENDDO
!
! saving the 2D data on the file if lwrite_loc=.t.
!
IF (lwrite_loc) THEN
IF(fil_loc.EQ.' ') CALL errore ('local','fil_loc no name',1)
CALL seqopn(4,fil_loc,'unformatted',exst)
WRITE(4) n2d
WRITE(4) ((newbg(ig,il), ig=1, ngper*npol), il=1, n2d)
#ifdef __PARA
CALL mp_barrier()
#endif
WRITE(4) (((psiper(ig,il,k),ig=1,n2d),il=1,n2d), &
k=1,nrzp)
WRITE(4) ((zkr(il,k),il=1,n2d),k=1,nrzp)
CLOSE(unit=4)
ENDIF
deallocate(amat)
deallocate(amat1)
deallocate(ymat)
deallocate(gp)
deallocate(fftxy)
DEALLOCATE(amat)
DEALLOCATE(amat1)
DEALLOCATE(ymat)
DEALLOCATE(gp)
DEALLOCATE(psibase)
DEALLOCATE(psiprob)
DEALLOCATE(el)
DEALLOCATE(fftxy)
CALL stop_clock('local')
RETURN
END SUBROUTINE local
!-----------------------------------
return
end subroutine local_2
FUNCTION number(gp, at, fftxy, nrx, nry)
!
@ -285,10 +398,9 @@ FUNCTION number(gp, at, fftxy, nrx, nry)
! and write on output its fft position.
!
IMPLICIT NONE
INTEGER :: nrx, nry, fftxy(-(nrx-1)/2:nrx/2, -(nry-1)/2:nry/2), &
INTEGER :: nrx, nry, fftxy(-nrx:nrx, -nrx:nry), &
number, n1, n2
REAL(kind=KIND(0.d0)) :: gp(2), at(3,3), x1, x2
REAL(kind=KIND(0.d0)), PARAMETER :: eps=1.d-4
x1=gp(1)*at(1,1)+gp(2)*at(2,1)
x2=gp(1)*at(1,2)+gp(2)*at(2,2)

5
PWCOND/local_set.f90
View File

@ -5,13 +5,13 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
subroutine local_set(n1, n2, n3, n4, n5, n6, n7)
subroutine local_set(n1, n2, n3, n4, n5, n6, n7, n8)
!
! To set up the number of all the orbitals to be 0 for
! local potential calculations
!
implicit none
integer :: n1, n2, n3, n4, n5, n6, n7
integer :: n1, n2, n3, n4, n5, n6, n7, n8
n1 = 0
n2 = 0
@ -20,6 +20,7 @@ subroutine local_set(n1, n2, n3, n4, n5, n6, n7)
n5 = 0
n6 = 0
n7 = 0
n8 = 0
return
end subroutine local_set

68
PWCOND/poten.f90
View File

@ -10,7 +10,7 @@
!
#include "f_defs.h"
!
SUBROUTINE poten
SUBROUTINE poten(vppot,nrz,z)
!
! This subroutine computes the 2D Fourier components of the
! local potential in each slab.
@ -18,29 +18,26 @@ SUBROUTINE poten
USE pwcom
USE noncollin_module, ONLY : noncolin, npol
USE cond
USE mp_global, ONLY : me_pool
USE pfft, ONLY : npp
USE mp, ONLY : mp_bcast
USE io_global, ONLY : ionode_id
IMPLICIT NONE
INTEGER :: &
i, j, ij, ijx, k, n, p, il, ik, kstart, klast, &
ix, jx, kx, ir, ir1, ixy, info
ix, jx, kx, ir, ir1, ixy, nrz, info
INTEGER :: iis, jjs, is(4), js(4), ispin, nspin_eff
INTEGER :: ionode_id
INTEGER, ALLOCATABLE :: ipiv(:)
REAL(kind=DP), PARAMETER :: eps = 1.d-8
REAL(kind=DP) :: arg, bet
REAL(kind=DP), ALLOCATABLE :: gz(:), auxr(:)
REAL(kind=DP) :: arg, bet, z(nrz+1), zlen
REAL(kind=DP), ALLOCATABLE :: gz(:), allv(:), auxr(:)
COMPLEX(kind=DP), PARAMETER :: cim = (0.d0,1.d0)
COMPLEX(kind=DP) :: caux
COMPLEX(kind=DP) :: caux, vppot(nrz,nrx*nry,npol,npol)
COMPLEX(kind=DP), ALLOCATABLE :: aux(:), amat(:,:), amat0(:,:)
COMPLEX(kind=DP), ALLOCATABLE :: vppot0(:,:,:,:)
LOGICAL :: lg
CALL start_clock('poten')
ALLOCATE( ipiv( nrz ) )
ALLOCATE( gz( nrz ) )
@ -49,6 +46,9 @@ SUBROUTINE poten
ALLOCATE( amat( nrz, nrz ) )
ALLOCATE( amat0( nrz, nrz ) )
zlen = at(3,3)
!
! Compute the Gz vectors in the z direction
!
@ -66,9 +66,9 @@ DO n=1,nrz
bet=gz(n)*(z(p+1)-z(p))*tpi
IF (ABS(gz(n)).GT.eps) THEN
caux=cim*(CMPLX(COS(bet),-SIN(bet))-(1.d0,0.d0)) &
/zl/gz(n)/tpi
/zlen/gz(n)/tpi
ELSE
caux=(z(p+1)-z(p))/zl
caux=(z(p+1)-z(p))/zlen
ENDIF
amat0(n,p)=CMPLX(COS(arg),-SIN(arg))*caux
ENDDO
@ -91,6 +91,11 @@ ENDIF
!
! To form local potential on the real space mesh
!
!
#ifdef __PARA
allocate ( allv(nrx1*nrx2*nrx3) )
#endif
vppot = 0.d0
DO ispin=1,nspin_eff
IF (noncolin) THEN
@ -105,37 +110,13 @@ DO ispin=1,nspin_eff
!
! To collect the potential from different CPUs
!
aux(:) = (0.d0,0.d0)
#ifdef __PARA
kstart = 1
DO i=1, me_pool
kstart = kstart+npp(i)
ENDDO
klast = kstart+npp(me_pool+1)-1
call gather( auxr, allv )
CALL mp_bcast( allv, ionode_id )
aux = CMPLX(allv)
#else
aux = CMPLX(auxr)
#endif
DO i = 1, nrx1*nrx2
DO k = 1, nr3
lg = .TRUE.
ir = i+(k-1)*nrx2*nrx1
ir1 = ir
#ifdef __PARA
IF(k.GE.kstart.AND.k.LE.klast) THEN
lg = .TRUE.
ir1 = i+(k-kstart)*nrx2*nrx1
ELSE
lg = .FALSE.
ENDIF
#endif
IF (lg) THEN
aux(ir) = auxr(ir1)
ENDIF
ENDDO
ENDDO
#ifdef __PARA
CALL reduce (2*nrx1*nrx2*nrx3,aux)
#endif
!
! To find FFT of the local potential
!
@ -199,7 +180,7 @@ IF (noncolin) THEN
ENDIF
! do p = 1, nrz
! write(6,'(i5,2f12.6)') p, real(vppot(p,105,2,2)), imag(vppot(p,105,2,2))
! write(6,'(i5,2f12.6)') p, real(vppot(p,1,1,1)), imag(vppot(p,1,1,1))
! enddo
! stop
@ -209,6 +190,9 @@ ENDIF
DEALLOCATE(auxr)
DEALLOCATE(amat)
DEALLOCATE(amat0)
#ifdef __PARA
deallocate(allv)
#endif
CALL stop_clock('poten')

40
PWCOND/rotproc.f90
View File

@ -11,7 +11,7 @@
#include "f_defs.h"
!
SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
fundl1, intw1, intw2, n2d, norbf, norbnow)
fundl1, intw1, intw2, n2d, norbf, norb)
!
! This subroutine implements a matching procedure to construct
! local and nonlocal functions on the whole region from those computed
@ -30,15 +30,17 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
! So in this case there are 3 matching steps.
!
USE kinds, ONLY : DP
USE mp_global, ONLY: nproc, me_pool
USE kinds, ONLY : DP
USE noncollin_module, ONLY : npol
USE parallel_include
USE mp_global, ONLY : nproc, me_pool
IMPLICIT NONE
INTEGER :: k, ig, n, lam, lam1, iorb, iorb1, norbf, norbnow, n2d, &
INTEGER :: k, ig, n, lam, lam1, iorb, iorb1, norbf, norb, n2d, &
ibound, numb, ninsl, ib, icolor, ikey, new_comm, info
INTEGER, ALLOCATABLE :: ipiv(:)
COMPLEX(kind=DP) :: fun0(n2d, 2*n2d), & ! phi_n(0)
@ -63,8 +65,8 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
ALLOCATE( y( n2d ) )
ALLOCATE( amat( 2*n2d, 2*n2d ) )
ALLOCATE( amat_aux( 2*n2d, 2*n2d ) )
ALLOCATE( vec( 2*n2d, 2*n2d+npol*norbnow ) )
ALLOCATE( vec_aux( 2*n2d, 2*n2d+npol*norbnow ) )
ALLOCATE( vec( 2*n2d, 2*n2d+npol*norb ) )
ALLOCATE( vec_aux( 2*n2d, 2*n2d+npol*norb ) )
ALLOCATE( ipiv( 2*n2d ) )
numb=0
@ -92,7 +94,7 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
vec(lam, lam1)=fun0(lam, lam1)
vec(n2d+lam, lam1)=fund0(lam, lam1)
ENDDO
DO iorb=1, npol*norbnow
DO iorb=1, npol*norb
vec(lam, 2*n2d+iorb)=funl0(lam, iorb)
vec(n2d+lam, 2*n2d+iorb)=fundl0(lam, iorb)
ENDDO
@ -107,7 +109,7 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
vec(lam, n2d+lam1)=-fun1(lam, n2d+lam1)
vec(n2d+lam, n2d+lam1)=-fund1(lam, n2d+lam1)
ENDDO
DO iorb=1, npol*norbnow
DO iorb=1, npol*norb
vec(lam, 2*n2d+iorb)=-funl1(lam, iorb)
vec(n2d+lam, 2*n2d+iorb)=-fundl1(lam, iorb)
ENDDO
@ -116,11 +118,11 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
ENDIF
CALL mpi_allreduce(amat, amat_aux, 2*2*n2d*2*n2d, MPI_REAL8, &
MPI_SUM, new_comm, info)
CALL mpi_allreduce(vec, vec_aux, 2*2*n2d*(2*n2d+npol*norbnow), &
CALL mpi_allreduce(vec, vec_aux, 2*2*n2d*(2*n2d+npol*norb), &
MPI_REAL8, MPI_SUM, new_comm, info)
CALL DCOPY(2*2*n2d*2*n2d, amat_aux, 1, amat, 1)
CALL DCOPY(2*2*n2d*(2*n2d+npol*norbnow), vec_aux, 1, vec, 1)
CALL ZGESV(2*n2d, 2*n2d+npol*norbnow, amat, 2*n2d, ipiv, &
CALL DCOPY(2*2*n2d*(2*n2d+npol*norb), vec_aux, 1, vec, 1)
CALL ZGESV(2*n2d, 2*n2d+npol*norb, amat, 2*n2d, ipiv, &
vec, 2*n2d, info)
!
! recalculate the functions for CPU which is left to matching
@ -136,7 +138,7 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
vec(n2d+lam, n)*fund1(ig, n2d+lam)
ENDDO
ENDDO
DO iorb=1, npol*norbnow
DO iorb=1, npol*norb
DO lam=1, n2d
funl1(ig, iorb)=funl1(ig, iorb)+ &
vec(n2d+lam, 2*n2d+iorb)*fun1(ig, n2d+lam)
@ -174,7 +176,7 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
vec(lam, n2d+n)*fund0(ig, lam)
ENDDO
ENDDO
DO iorb=1, npol*norbnow
DO iorb=1, npol*norb
DO lam=1, n2d
funl0(ig, iorb)=funl0(ig, iorb)+ &
vec(lam, 2*n2d+iorb)*fun0(ig, lam)
@ -202,15 +204,15 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
! to recalculate the integrals for a given group of CPU
!
IF((me_pool+1).GE.ib) THEN
DO iorb=1, npol*norbnow
DO iorb1=1, npol*norbnow
DO iorb=1, npol*norb
DO iorb1=1, npol*norb
DO lam=1, n2d
intw2(iorb,iorb1)=intw2(iorb,iorb1)+ &
vec(n2d+lam, 2*n2d+iorb1)*intw1(iorb, n2d+lam)
ENDDO
ENDDO
ENDDO
DO iorb=1, npol*norbnow
DO iorb=1, npol*norb
x=(0.d0,0.d0)
DO n=1, n2d
DO lam=1, n2d
@ -224,15 +226,15 @@ SUBROUTINE rotproc (fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
ENDDO
ENDDO
ELSE
DO iorb=1, npol*norbnow
DO iorb1=1, npol*norbnow
DO iorb=1, npol*norb
DO iorb1=1, npol*norb
DO lam=1, n2d
intw2(iorb, iorb1)=intw2(iorb, iorb1)+ &
vec(lam, 2*n2d+iorb1)*intw1(iorb, lam)
ENDDO
ENDDO
ENDDO
DO iorb=1, npol*norbnow
DO iorb=1, npol*norb
x=(0.d0,0.d0)
DO n=1, n2d
DO lam=1, n2d

345
PWCOND/save_cond.f90 Normal file
View File

@ -0,0 +1,345 @@
!
! 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 save_cond (lwrite, lsr, ef, nrz, nocros, noins, &
norb, r, rab, betar)
!
! This subroutine writes/reads variables needed for PWCOND
! so that the punch file from PW calculations is not needed.
!
use kinds, only : DP
USE parameters, only : ndmx, nbrx, npsx
USE cell_base, ONLY : alat, tpiba, tpiba2, at, bg
use io_files, only : save_file
use lsda_mod, only: nspin
USE noncollin_module, ONLY : noncolin, npol
use spin_orb, only : lspinorb
use cond, only : sarea, nrx, nry, norbf, tblml, crosl, taunewl, &
zpseul, zpseul_nc, zl, vppotl, tblms, cross, taunews, zpseus,&
zpseus_nc, zs, vppots, tblmr, crosr, taunewr, zpseur, &
zpseur_nc, zr, vppotr, iofspin
implicit none
integer :: lsr, nrz, nocros, noins, norb, i, j, k, l, m
logical :: lwrite
REAL(kind=DP) :: ef, r(ndmx,npsx), rab(ndmx,npsx), &
betar(ndmx,nbrx,npsx)
integer, ALLOCATABLE :: ind(:,:), tblm(:,:), cros(:,:)
REAL(kind=DP), ALLOCATABLE :: z(:), zpseu(:,:,:), re(:,:,:,:), &
im(:,:,:,:), c(:), taunew(:,:)
COMPLEX(kind=DP), ALLOCATABLE :: vppot(:,:,:,:), zpseu_nc(:,:,:,:)
character(len=2) :: ext
call start_clock('save_cond')
if(lsr.eq.1) then
ext='.l'
elseif(lsr.eq.2) then
ext='.s'
elseif(lsr.eq.3) then
ext='.r'
endif
if (lwrite) then
allocate( vppot(nrz, nrx * nry, npol, npol) )
allocate( z(nrz+1) )
allocate( taunew(4,norb) )
allocate( tblm(4,norb) )
allocate( cros(norb, nrz) )
allocate( zpseu(2, norb, norb) )
if (noncolin) allocate(zpseu_nc(2, norb, norb, nspin))
if(lsr.eq.1) then
vppot = vppotl
z = zl
taunew = taunewl
tblm = tblml
cros = crosl
zpseu = zpseul
if (noncolin) zpseu_nc = zpseul_nc
elseif(lsr.eq.2) then
vppot = vppots
z = zs
taunew = taunews
tblm = tblms
cros = cross
zpseu = zpseus
if (noncolin) zpseu_nc = zpseus_nc
elseif(lsr.eq.3) then
vppot = vppotr
z = zr
taunew = taunewr
tblm = tblmr
cros = crosr
zpseu = zpseur
if (noncolin) zpseu_nc = zpseur_nc
endif
open (3,file=trim(save_file)//ext,form='formatted', &
status='unknown')
write(3,*) nspin, npol, noncolin, lspinorb
if(nspin.eq.2) write(3,*) iofspin
write(3,*) alat, tpiba, tpiba2
write(3,'(6f20.14)') ((at(i,j),i=1,3),j=1,3)
write(3,'(6f20.14)') ((bg(i,j),i=1,3),j=1,3)
write(3,*) sarea
write(3,*) ef
write(3,*) nrx, nry, nrz
write(3,*) nocros, noins, norb, norbf
write(3,'(40i3)') ((tblm(i,j),i=1,4),j=1,norb)
write(3,'(120i1)') ((cros(j,i),i=1,nrz),j=1,norb)
write(3,'(6f20.14)') ((taunew(i,j),i=1,4),j=1,norb)
! write zpseu
allocate( ind(3,2*norb*norb) )
allocate( c(2*norb*norb) )
m=0
do i=1, 2
do j=1, norb
do k=1, norb
if(abs(zpseu(i,j,k)).gt.1.d-12) then
m = m+1
ind(1,m) = i
ind(2,m) = j
ind(3,m) = k
c(m) = zpseu(i,j,k)
endif
enddo
enddo
enddo
write(3,*) m
write(3,'(25i5)') ((ind(i,j),i=1,3),j=1,m)
write(3,'(6f20.14)') (c(i),i=1,m)
deallocate(ind)
deallocate(c)
if(noncolin) then
write(3,'(6f20.14)') ((((real(zpseu_nc(i,j,k,l)),i=1,2), &
j=1,norb),k=1,norb),l=1,nspin)
write(3,'(6f20.14)') ((((aimag(zpseu_nc(i,j,k,l)),i=1,2), &
j=1,norb),k=1,norb),l=1,nspin)
endif
write(3,'(6f20.14)') (z(i), i=1, nrz+1)
write(3,'(6f20.14)') ((((real(vppot(i,j,k,l)),i=1,nrz), &
j=1,nrx*nry),k=1,npol),l=1,npol)
write(3,'(6f20.14)') ((((aimag(vppot(i,j,k,l)),i=1,nrz), &
j=1,nrx*nry),k=1,npol),l=1,npol)
! write r
allocate( ind(2,npsx*ndmx) )
allocate( c(npsx*ndmx) )
m=0
do i=1, ndmx
do j=1, npsx
if(abs(r(i,j)).gt.1.d-12) then
m = m+1
ind(1,m) = i
ind(2,m) = j
c(m) = r(i,j)
endif
enddo
enddo
write(3,*) m
write(3,'(25i5)') ((ind(i,j),i=1,2),j=1,m)
write(3,'(6f20.14)') (c(i),i=1,m)
deallocate(ind)
deallocate(c)
! write rab
allocate( ind(2,npsx*ndmx) )
allocate( c(npsx*ndmx) )
m=0
do i=1, ndmx
do j=1, npsx
if(abs(rab(i,j)).gt.1.d-12) then
m = m+1
ind(1,m) = i
ind(2,m) = j
c(m) = rab(i,j)
endif
enddo
enddo
write(3,*) m
write(3,'(25i5)') ((ind(i,j),i=1,2),j=1,m)
write(3,'(6f20.14)') (c(i),i=1,m)
deallocate(ind)
deallocate(c)
! write betar
allocate( ind(3,npsx*nbrx*ndmx) )
allocate( c(npsx*nbrx*ndmx) )
m=0
do i=1, ndmx
do j=1, nbrx
do k=1, npsx
if(abs(betar(i,j,k)).gt.1.d-12) then
m = m+1
ind(1,m) = i
ind(2,m) = j
ind(3,m) = k
c(m) = betar(i,j,k)
endif
enddo
enddo
enddo
write(3,*) m
write(3,'(25i5)') ((ind(i,j),i=1,3),j=1,m)
write(3,'(6f20.14)') (c(i),i=1,m)
deallocate(ind)
deallocate(c)
close(unit=3)
else
open (3,file=trim(save_file)//ext,form='formatted', &
status='unknown')
read(3,*) nspin, npol, noncolin, lspinorb
if(nspin.eq.2) read(3,*) iofspin
read(3,*) alat, tpiba, tpiba2
read(3,'(6f20.14)') ((at(i,j),i=1,3),j=1,3)
read(3,'(6f20.14)') ((bg(i,j),i=1,3),j=1,3)
read(3,*) sarea
read(3,*) ef
read(3,*) nrx, nry, nrz
read(3,*) nocros, noins, norb, norbf
allocate( vppot(nrz, nrx * nry, npol, npol) )
allocate( z(nrz+1) )
allocate( taunew(4,norb) )
allocate( tblm(4,norb) )
allocate( cros(norb, nrz) )
allocate( zpseu(2, norb, norb) )
if (noncolin) allocate(zpseu_nc(2, norb, norb, nspin))
read(3,'(40i3)') ((tblm(i,j),i=1,4),j=1,norb)
read(3,'(120i1)') ((cros(j,i),i=1,nrz),j=1,norb)
read(3,'(6f20.14)') ((taunew(i,j),i=1,4),j=1,norb)
! read zpseu
read(3,*) m
allocate( ind(3,m) )
allocate( c(m) )
read(3,'(25i5)') ((ind(i,j),i=1,3),j=1,m)
read(3,'(6f20.14)') (c(i),i=1,m)
zpseu = 0.d0
do i=1, m
zpseu(ind(1,i),ind(2,i),ind(3,i)) = c(i)
enddo
deallocate(ind)
deallocate(c)
if(noncolin) then
allocate ( re(2,norb,norb,nspin) )
allocate ( im(2,norb,norb,nspin) )
read(3,'(6f20.14)') ((((re(i,j,k,l),i=1,2), &
j=1,norb),k=1,norb),l=1,nspin)
read(3,'(6f20.14)') ((((im(i,j,k,l),i=1,2), &
j=1,norb),k=1,norb),l=1,nspin)
zpseu_nc = CMPLX(re,im)
deallocate(re)
deallocate(im)
endif
read(3,'(6f20.14)') (z(i), i=1, nrz+1)
allocate ( re(nrz,nrx*nry,npol,npol) )
allocate ( im(nrz,nrx*nry,npol,npol) )
read(3,'(6f20.14)') ((((re(i,j,k,l),i=1,nrz), &
j=1,nrx*nry),k=1,npol),l=1,npol)
read(3,'(6f20.14)') ((((im(i,j,k,l),i=1,nrz), &
j=1,nrx*nry),k=1,npol),l=1,npol)
vppot = CMPLX(re,im)
deallocate(re)
deallocate(im)
! read r
read(3,*) m
allocate( ind(2,m) )
allocate( c(m) )
read(3,'(25i5)') ((ind(i,j),i=1,2),j=1,m)
read(3,'(6f20.14)') (c(i),i=1,m)
r = 0.d0
do i=1,m
r(ind(1,i),ind(2,i))=c(i)
enddo
deallocate(ind)
deallocate(c)
! read rab
read(3,*) m
allocate( ind(2,m) )
allocate( c(m) )
read(3,'(25i5)') ((ind(i,j),i=1,2),j=1,m)
read(3,'(6f20.14)') (c(i),i=1,m)
rab = 0.d0
do i=1,m
rab(ind(1,i),ind(2,i))=c(i)
enddo
deallocate(ind)
deallocate(c)
! read betar
read(3,*) m
allocate( ind(3,m) )
allocate( c(m) )
read(3,'(25i5)') ((ind(i,j),i=1,3),j=1,m)
read(3,'(6f20.14)') (c(i),i=1,m)
betar = 0.d0
do i=1,m
betar(ind(1,i),ind(2,i),ind(3,i))=c(i)
enddo
deallocate(ind)
deallocate(c)
close(unit=3)
if(lsr.eq.1) then
allocate( vppotl(nrz, nrx * nry, npol, npol) )
allocate( zl(nrz+1) )
allocate( taunewl(4,norb) )
allocate( tblml(4,norb) )
allocate( crosl(norb, nrz) )
allocate( zpseul(2, norb, norb) )
if (noncolin) allocate(zpseul_nc(2, norb, norb, nspin))
vppotl = vppot
zl = z
taunewl = taunew
tblml = tblm
crosl = cros
zpseul = zpseu
if (noncolin) zpseul_nc = zpseu_nc
elseif(lsr.eq.2) then
allocate( vppots(nrz, nrx * nry, npol, npol) )
allocate( zs(nrz+1) )
allocate( taunews(4,norb) )
allocate( tblms(4,norb) )
allocate( cross(norb, nrz) )
allocate( zpseus(2, norb, norb) )
if (noncolin) allocate(zpseus_nc(2, norb, norb, nspin))
vppots = vppot
zs = z
taunews = taunew
tblms = tblm
cross = cros
zpseus = zpseu
if (noncolin) zpseus_nc = zpseu_nc
elseif(lsr.eq.3) then
allocate( vppotr(nrz, nrx * nry, npol, npol) )
allocate( zr(nrz+1) )
allocate( taunewr(4,norb) )
allocate( tblmr(4,norb) )
allocate( crosr(norb, nrz) )
allocate( zpseur(2, norb, norb) )
if (noncolin) allocate(zpseur_nc(2, norb, norb, nspin))
vppotr = vppotr
zr = z
taunewr = taunew
tblmr = tblm
crosr = cros
zpseur = zpseu
if (noncolin) zpseur_nc = zpseu_nc
endif
endif
deallocate( vppot )
deallocate( z )
deallocate( taunew )
deallocate( tblm )
deallocate( cros )
deallocate( zpseu )
if (noncolin) deallocate(zpseu_nc)
call stop_clock('save_cond')
end subroutine save_cond

49
PWCOND/scatter_back.f90
View File

@ -1,4 +1,3 @@
!
! Copyright (C) 2003 A. Smogunov
! This file is distributed under the terms of the
@ -9,27 +8,27 @@
! Generalized to spinor wavefunctions and spin-orbit Oct. 2004 (ADC).
!
!
subroutine scatter_back(app, bpp, an, bn, af, ci, di, ezk, emzk, &
s1m, s2m, s3m, s4m, &
startk, lastk, norbnow, orbin, dz)
subroutine scatter_back(psiper, zk, app, bpp, an, bn, af, ci, di, &
ezk, emzk, s1m, s2m, s3m, s4m, kin, kfin, &
nrz, nrzp, norb, cros, dz)
!
! This subroutine computes the second half of the local functions
! and their integrals.
!
#include "f_defs.h"
use pwcom
use noncollin_module, ONLY : npol
use cond
implicit none
integer :: kz, n, lam, lam1, k, iorb, startk, lastk, norbnow, &
orbin, iorba
integer :: kp, n, lam, lam1, k, iorb, kin, kfin, norb, &
nrz, nrzp, iorba, cros(norb,nrz)
real(kind=DP) :: dz
complex(kind=DP), parameter :: cim=(0.d0, 1.d0), one=(1.d0,0.d0)
complex(kind=DP) :: c, s1, s2, s3, s4, ZDOTC, CONJG, &
psiper(n2d,n2d,nrzp), zk(n2d,nrzp), &
app(n2d,n2d), an(n2d,n2d), bn(n2d,n2d), &
bpp(n2d,n2d), af(n2d,n2d), &
ci(norbnow*npol,n2d,nrzp), di(norbnow*npol,n2d,nrzp), &
ci(norb*npol,n2d,nrzp), di(norb*npol,n2d,nrzp), &
ezk(n2d,nrzp), emzk(n2d,nrzp), &
s1m(n2d,n2d),s2m(n2d,n2d),s3m(n2d,n2d),s4m(n2d,n2d)
@ -46,13 +45,12 @@ subroutine scatter_back(app, bpp, an, bn, af, ci, di, ezk, emzk, &
af(lam,lam)=(1.d0,0.d0)
app(lam,lam)=(1.d0,0.d0)
enddo
kz=nkofz(lastk)
do iorb=1, norbnow*npol
do iorb=1, norb*npol
iorba=iorb
if (npol.eq.2) iorba=(iorb+1)/2
if (cross(orbin+iorba,lastk).eq.1) then
if (cros(iorba,kfin).eq.1) then
do n=1, n2d
intw1(iorb,n2d+n)=ci(iorb,n,kz)
intw1(iorb,n2d+n)=ci(iorb,n,nrzp)
enddo
endif
enddo
@ -60,14 +58,14 @@ subroutine scatter_back(app, bpp, an, bn, af, ci, di, ezk, emzk, &
!
! The loop over slabs
!
do k=lastk-1, startk, -1
kz=nkofz(k)
do k=kfin-1, kin, -1
kp = k-kin+1
do lam=1, n2d
do lam1=1,n2d
c=ZDOTC(n2d,psiper(1,lam,kz),1,psiper(1,lam1,kz+1),1)
s1m(lam,lam1)=(zk(lam,kz)+zk(lam1,kz+1))/zk(lam,kz)*c
s2m(lam,lam1)=(zk(lam,kz)-zk(lam1,kz+1))/zk(lam,kz)*c
c=ezk(lam1,kz+1)
c=ZDOTC(n2d,psiper(1,lam,kp),1,psiper(1,lam1,kp+1),1)
s1m(lam,lam1)=(zk(lam,kp)+zk(lam1,kp+1))/zk(lam,kp)*c
s2m(lam,lam1)=(zk(lam,kp)-zk(lam1,kp+1))/zk(lam,kp)*c
c=ezk(lam1,kp+1)
s3m(lam,lam1)=s1m(lam,lam1)*c
s4m(lam,lam1)=s2m(lam,lam1)*c
enddo
@ -80,18 +78,18 @@ subroutine scatter_back(app, bpp, an, bn, af, ci, di, ezk, emzk, &
bn=bn*0.5d0
do lam=1,n2d
do n=1, n2d
an(lam,n)=an(lam,n)*emzk(lam,kz)
an(lam,n)=an(lam,n)*emzk(lam,kp)
enddo
enddo
do iorb=1, norbnow*npol
do iorb=1, norb*npol
iorba=iorb
if (npol.eq.2) iorba=(iorb+1)/2
if (cross(orbin+iorba,k).eq.1) then
if (cros(iorba,k).eq.1) then
do lam=1, n2d
do n=1, n2d
intw1(iorb,n2d+n)=intw1(iorb,n2d+n)+ &
an(lam,n)*ci(iorb,lam,kz)+ &
bn(lam,n)*di(iorb,lam,kz)
an(lam,n)*ci(iorb,lam,kp)+ &
bn(lam,n)*di(iorb,lam,kp)
enddo
enddo
endif
@ -100,7 +98,10 @@ subroutine scatter_back(app, bpp, an, bn, af, ci, di, ezk, emzk, &
call DCOPY(2*n2d*n2d, bn, 1, bpp, 1)
an=(0.d0,0.d0)
bn=(0.d0,0.d0)
call rotateb(app, bpp, af, intw1, n2d, norbf, norbnow, npol)
call rotateb(app, bpp, af, intw1, n2d, norbf, norb, npol)
! write(6,*) k, kin
enddo
call stop_clock('scatter_back')

383
PWCOND/scatter_forw.f90
View File

@ -11,7 +11,8 @@
!
#include "f_defs.h"
!
SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
subroutine scatter_forw(nrz, nrzp, z, psiper, zk, norb, tblm, cros, &
taunew, r, rab, betar)
!
! This subroutine computes local Phi_n and partial nonlocal Phi_alpha
! solutions of the Schrodinger equation in the region zin<z<zfin
@ -20,31 +21,31 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
! It computes also the integrals (intw1, intw2) of Phi_n and
! Phi_alpha over beta-functions inside the unit cell.
!
USE pwcom
USE constants, ONLY : tpi
USE parameters, only : ndmx, nbrx, npsx
USE cell_base, ONLY : tpiba
USE noncollin_module, ONLY : npol
USE cond
!
IMPLICIT NONE
INTEGER :: &
orbin, & ! starting orbital in zin<z<zfin
orbfin, & ! final orbital in zin<z<zfin
norbnow, & ! total number of orbitals in zin<z<zfin
kin, & ! first slab in zin<z<zfin
kfin, & ! last slab in zin<z<zfin
startk, & ! first slab for a given CPU
lastk, & ! last slab for a given CPU
nrz, nrzp, norb,&
cros(norb,nrz), &
tblm(4,norb), &
k, kz, n, lam, ig, lam1, mdim, itt, nbb, iorb, iorb1, &
iorbs, iorb1s, iorba, iorb1a, is, kkz, nok
iorba, iorb1a, is, kp, nok, k1, nt, nb, kin, kfin
INTEGER :: info
INTEGER, ALLOCATABLE :: inslab(:)
real(kind=DP) :: z(nrz+1), r(1:ndmx,npsx), rab(1:ndmx,npsx), &
betar(1:ndmx,nbrx,npsx), taunew(4,norb)
REAL(kind=DP), PARAMETER :: eps=1.d-8
REAL(kind=DP) :: DDOT, zin, zfin, dz, tr, tr1, dz1
REAL(kind=DP) :: DDOT, dz, tr, tr1, dz1
COMPLEX(kind=DP), PARAMETER :: cim=(0.d0,1.d0), one=(1.d0, 0.d0), &
zero=(0.d0,0.d0)
COMPLEX(kind=DP) :: int1d, int2d, c, d, e, f, s1, s2, s3, s4, arg,&
f1p, ZDOTC, fact, factm
f1p, ZDOTC, fact, factm, psiper(n2d,n2d,nrzp), &
zk(n2d,nrzp)
COMPLEX(kind=DP), ALLOCATABLE :: &
psigper(:,:), & ! psigper(g,lam)=newbg(g,lam1) psiper(lam1,lam)
w0(:,:,:), & ! w0(z,g,m) are 2D Fourier components (see four.f)
@ -66,89 +67,80 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
ezk1(:,:), emzk1(:,:)
CALL start_clock('scatter_forw')
norbnow = orbfin - orbin + 1
orbin = orbin - 1
!
! Find first and last slab (kin, kfin)
!
DO k=1, nrz
IF (z(k).LE.zin+eps) kin=k
IF (z(k).LE.zfin-eps) kfin=k
ENDDO
dz=z(kin+1)-z(kin)
dz1=dz/nz1
dz = z(2) - z(1)
dz1 = dz/nz1
!
! Divide the slabs among CPU
!
CALL divide(kfin-kin+1,startk,lastk)
startk=kin+startk-1
lastk=kin+lastk-1
call divide(nrz,kin,kfin)
!------------------------
! Start of 2D Fourier components calculations and depending
! variables
!
ALLOCATE( psigper( ngper*npol, n2d ) )
ALLOCATE( w( nz1, n2d, norbnow*npol ) )
ALLOCATE( w( nz1, n2d, norb*npol ) )
ALLOCATE( w0( nz1, ngper, 5 ) )
ALLOCATE( cix( nz1, n2d, norbnow*npol ) )
ALLOCATE( dix( nz1, n2d, norbnow*npol ) )
ALLOCATE( ci( norbnow*npol, n2d, nrzp ) )
ALLOCATE( di( norbnow*npol, n2d, nrzp ) )
ALLOCATE( inslab( norbnow ) )
ALLOCATE( cix( nz1, n2d, norb*npol ) )
ALLOCATE( dix( nz1, n2d, norb*npol ) )
ALLOCATE( ci( norb*npol, n2d, nrzp ) )
ALLOCATE( di( norb*npol, n2d, nrzp ) )
ALLOCATE( inslab( norb ) )
ALLOCATE( ezk( n2d, nrzp ) )
ALLOCATE( emzk( n2d, nrzp ) )
ALLOCATE( ezk1( nz1, n2d ) )
ALLOCATE( emzk1( nz1, n2d ) )
ALLOCATE( zk2( n2d, nrzp ) )
intw1=(0.d0,0.d0)
intw2=(0.d0,0.d0)
DO k=startk,lastk
kz=nkofz(k)
DO lam=1,n2d
arg=cim*tpi*zk(lam, kz)*dz
ezk(lam,kz)=EXP(arg)
emzk(lam,kz)=EXP(-arg)
arg=cim*tpi*zk(lam, kz)*dz1
zk2(lam,kz)=cim/(2.d0*zk(lam,kz)*tpiba)
ENDDO
ENDDO
do k=kin,kfin
kp = k-kin+1
do lam=1,n2d
arg=cim*tpi*zk(lam, kp)*dz
ezk(lam,kp)=exp(arg)
emzk(lam,kp)=exp(-arg)
zk2(lam,kp)=cim/(2.d0*zk(lam,kp)*tpiba)
enddo
enddo
!
! some orbitals relations
!
DO iorb=1, norbnow
inslab(iorb)=0
ENDDO
DO iorb=1, norbnow
iorbs=orbin+iorb
IF(inslab(iorb).EQ.0.AND.mnew(iorbs).EQ.1) THEN
itt=itnew(iorbs)
nbb=nbnew(iorbs)
DO iorb1=iorb, norbnow
iorb1s=orbin+iorb1
IF(mnew(iorb1s).EQ.1.AND.nbnew(iorb1s).EQ.nbb) THEN
tr=ABS(taunew(3,iorbs)-taunew(3,iorb1s))
IF(itnew(iorb1s).EQ.itt.AND.tr.LE.eps) THEN
inslab(iorb1)=iorb
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
do iorb=1, norb
inslab(iorb) = 0
enddo
do iorb = 1, norb
if(inslab(iorb).eq.0.and.tblm(4,iorb).eq.1) then
itt = tblm(1,iorb)
nbb = tblm(2,iorb)
do iorb1 = iorb, norb
if(tblm(4,iorb1).eq.1.and.tblm(2,iorb1).eq.nbb) then
tr = abs(taunew(3,iorb)-taunew(3,iorb1))
if(tblm(1,iorb1).eq.itt.and.tr.le.eps) then
inslab(iorb1) = iorb
endif
endif
enddo
endif
enddo
CALL start_clock('integrals')
!
! The loop over slabs to compute ci, di, and initial intw2
!
DO k=startk, lastk
do k = kin, kfin
kp = k-kin+1
! write(6,*) 'integrals k=', k
kkz=nkofz(k)
DO lam=1,n2d
arg=cim*zk(lam,kkz)*dz1*tpi
arg=cim*zk(lam,kp)*dz1*tpi
fact=EXP(arg)
factm=EXP(-arg)
ezk1(1,lam)=fact
@ -159,19 +151,26 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
ENDDO
ENDDO
CALL ZGEMM('n', 'n', ngper*npol, n2d, n2d, one, newbg, ngper*npol, &
psiper(1,1,kkz), n2d, zero, psigper, ngper*npol)
if(ewind.le.100.d0) then
CALL ZGEMM('n', 'n', ngper*npol, n2d, n2d, one, newbg, &
ngper*npol, psiper(1,1,kp), n2d, zero, &
psigper, ngper*npol)
else
psigper(:,:) = psiper(:,:,kp)
endif
w=(0.d0,0.d0)
DO iorb=1, norbnow
iorbs=orbin+iorb
IF(cross(iorbs,k).EQ.1.AND.inslab(iorb).EQ.iorb) THEN
mdim=2*ls(iorbs)+1
CALL four(iorbs, w0, k, dz)
DO iorb1=1, norbnow
iorb1s=orbin+iorb1
DO iorb=1, norb
IF(cros(iorb,k).EQ.1.AND.inslab(iorb).EQ.iorb) THEN
mdim = 2*tblm(3,iorb)+1
nt = tblm(1,iorb)
nb = tblm(2,iorb)
call four(w0, z(k), dz, tblm(1,iorb), taunew(1,iorb), &
r(1,nt), rab(1,nt), betar(1,nb,nt))
DO iorb1=1, norb
IF(inslab(iorb1).EQ.iorb) THEN
DO ig=1, ngper
tr=-tpi*DDOT(2,gper(1,ig),1,taunew(1,iorb1s),1)
tr=-tpi*DDOT(2,gper(1,ig),1,taunew(1,iorb1),1)
c=CMPLX(COS(tr),SIN(tr))
DO lam=1, n2d
DO is=1,npol
@ -189,38 +188,35 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
ENDDO
ENDIF
ENDDO
DO iorb=1, norbnow*npol
DO iorb=1, norb*npol
iorba=iorb
IF (npol.EQ.2) iorba=(iorb+1)/2
iorbs = orbin + iorba
IF (cross(iorbs,k).EQ.1) THEN
IF (cros(iorba,k).EQ.1) THEN
DO lam=1, n2d
CALL setint(w(1,lam,iorb),cix(1,lam,iorb),dix(1,lam,iorb), &
ezk1(1,lam), emzk1(1,lam), nz1)
ENDDO
ENDIF
ENDDO
DO iorb=1, norbnow*npol
DO iorb=1, norb*npol
iorba=iorb
IF (npol.EQ.2) iorba=(iorb+1)/2
iorbs = orbin + iorba
IF (cross(iorbs,k).EQ.1) THEN
IF (cros(iorba,k).EQ.1) THEN
DO lam=1, n2d
ci(iorb,lam,kkz)=int1d(w(1,lam,iorb), &
zk(lam,kkz),dz,dz1,nz1,tpiba,1)
di(iorb,lam,kkz)=int1d(w(1,lam,iorb), &
zk(lam,kkz),dz,dz1,nz1,tpiba,-1)
ci(iorb,lam,kp)=int1d(w(1,lam,iorb), &
zk(lam,kp),dz,dz1,nz1,tpiba,1)
di(iorb,lam,kp)=int1d(w(1,lam,iorb), &
zk(lam,kp),dz,dz1,nz1,tpiba,-1)
ENDDO
DO iorb1=1, norbnow*npol
DO iorb1=1, norb*npol
iorb1a=iorb1
IF (npol.EQ.2) iorb1a=(iorb1+1)/2
iorb1s = orbin + iorb1a
IF (cross(iorb1s,k).EQ.1) THEN
IF (cros(iorb1a,k).EQ.1) THEN
DO lam=1, n2d
intw2(iorb,iorb1)=intw2(iorb,iorb1)- &
int2d(w(1,lam,iorb),w(1,lam,iorb1),cix(1,lam,iorb1), &
dix(1,lam,iorb1),ezk1(1,lam),emzk1(1,lam), &
zk(lam,kkz),dz1,tpiba,nz1)*zk2(lam,kkz)
zk(lam,kp),dz1,tpiba,nz1)*zk2(lam,kp)
ENDDO
ENDIF
ENDDO
@ -255,14 +251,15 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
ALLOCATE( s2m( n2d, n2d ) )
ALLOCATE( s3m( n2d, n2d ) )
ALLOCATE( s4m( n2d, n2d ) )
ALLOCATE( bnlf( n2d, norbnow*npol ) )
ALLOCATE( anln( n2d, norbnow*npol ) )
ALLOCATE( bnln( n2d, norbnow*npol ) )
ALLOCATE( anlp( n2d, norbnow*npol ) )
ALLOCATE( bnlp( n2d, norbnow*npol ) )
ALLOCATE( anll( n2d, norbnow*npol ) )
ALLOCATE( ff( n2d, norbnow*npol ) )
ALLOCATE( fl( n2d, norbnow*npol ) )
ALLOCATE( bnlf( n2d, norb*npol ) )
ALLOCATE( anln( n2d, norb*npol ) )
ALLOCATE( bnln( n2d, norb*npol ) )
ALLOCATE( anlp( n2d, norb*npol ) )
ALLOCATE( bnlp( n2d, norb*npol ) )
ALLOCATE( anll( n2d, norb*npol ) )
ALLOCATE( ff( n2d, norb*npol ) )
ALLOCATE( fl( n2d, norb*npol ) )
!
! We set up the starting values
@ -287,21 +284,19 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
!
! To compute intw1, ff, fl for the first slab
!
kz=nkofz(startk)
DO iorb=1, norbnow*npol
DO iorb=1, norb*npol
iorba=iorb
IF (npol.EQ.2) iorba=(iorb+1)/2
iorbs = orbin + iorba
IF (cross(iorbs, startk).EQ.1) THEN
IF (cros(iorba,kin).EQ.1) THEN
DO lam=1, n2d
intw1(iorb,lam)=di(iorb, lam, kz)
arg=ezk(lam, kz)
IF (ABS(DIMAG(zk(lam, kz))).LT.eps) THEN
ff(lam,iorb)=-arg*CONJG(di(iorb,lam,kz))*zk2(lam,kz)
fl(lam,iorb)=-arg*CONJG(ci(iorb,lam,kz))*zk2(lam,kz)
intw1(iorb,lam)=di(iorb, lam, 1)
arg=ezk(lam, 1)
IF (ABS(AIMAG(zk(lam, 1))).LT.eps) THEN
ff(lam,iorb)=-arg*CONJG(di(iorb,lam,1))*zk2(lam,1)
fl(lam,iorb)=-arg*CONJG(ci(iorb,lam,1))*zk2(lam,1)
ELSE
ff(lam,iorb)=-CONJG(ci(iorb,lam,kz))*zk2(lam,kz)
fl(lam,iorb)=-CONJG(di(iorb,lam,kz))*zk2(lam,kz)
ff(lam,iorb)=-CONJG(ci(iorb,lam,1))*zk2(lam,1)
fl(lam,iorb)=-CONJG(di(iorb,lam,1))*zk2(lam,1)
ENDIF
ENDDO
ENDIF
@ -310,26 +305,26 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
!------------------------------------
! The main loop over slabs
!
DO k=startk+1, lastk
CALL start_clock('scatter')
kz=nkofz(k)
DO iorb=1, norbnow*npol
DO k=kin+1, kfin
kp = k-kin+1
nok = 0
DO iorb=1, norb*npol
iorba=iorb
IF (npol.EQ.2) iorba=(iorb+1)/2
iorbs = orbin + iorba
tr=taunew(3,iorbs)-rsph(nbnew(iorbs),itnew(iorbs))
tr=taunew(3,iorba)-taunew(4,iorba)
IF (z(k)+dz.GT.tr) nok=iorb
ENDDO
DO lam=1, n2d
DO lam1=1,n2d
c=ZDOTC(n2d,psiper(1,lam,kz),1,psiper(1,lam1,kz-1),1)
s1m(lam,lam1)=(zk(lam,kz)+zk(lam1,kz-1))/zk(lam,kz)*c
s2m(lam,lam1)=(zk(lam,kz)-zk(lam1,kz-1))/zk(lam,kz)*c
c=ezk(lam1,kz-1)
c=ZDOTC(n2d,psiper(1,lam,kp),1,psiper(1,lam1,kp-1),1)
s1m(lam,lam1)=(zk(lam,kp)+zk(lam1,kp-1))/zk(lam,kp)*c
s2m(lam,lam1)=(zk(lam,kp)-zk(lam1,kp-1))/zk(lam,kp)*c
c=ezk(lam1,kp-1)
s3m(lam,lam1)=s1m(lam,lam1)*c
s4m(lam,lam1)=s2m(lam,lam1)*c
ENDDO
ENDDO
CALL ZGEMM('n','n',n2d,n2d,n2d,one,s3m,n2d,app,n2d,one,an,n2d)
CALL ZGEMM('n','n',n2d,n2d,n2d,one,s4m,n2d,app,n2d,one,bn,n2d)
an= an+s2m
@ -344,59 +339,57 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
bnln=bnln*0.5d0
DO lam=1, n2d
DO n=1, n2d
bn(lam,n)=bn(lam,n)*emzk(lam,kz)
bn(lam,n)=bn(lam,n)*emzk(lam,kp)
ENDDO
DO iorb=1, nok
bnln(lam,iorb)=bnln(lam,iorb)*emzk(lam,kz)
bnln(lam,iorb)=bnln(lam,iorb)*emzk(lam,kp)
ENDDO
ENDDO
CALL DCOPY(2*n2d*n2d, an, 1, app, 1)
CALL DCOPY(2*n2d*n2d, bn, 1, bpp, 1)
an=(0.d0,0.d0)
bn=(0.d0,0.d0)
CALL DCOPY(2*norbnow*npol*n2d, anln, 1, anlp, 1)
CALL DCOPY(2*norbnow*npol*n2d, bnln, 1, bnlp, 1)
CALL DCOPY(2*norb*npol*n2d, anln, 1, anlp, 1)
CALL DCOPY(2*norb*npol*n2d, bnln, 1, bnlp, 1)
anln=(0.d0,0.d0)
bnln=(0.d0,0.d0)
fl=(0.d0,0.d0)
!
DO iorb=1, norbnow*npol
DO iorb=1, norb*npol
iorba=iorb
IF (npol.EQ.2) iorba=(iorb+1)/2
iorbs = orbin + iorba
IF(cross(iorbs, k).EQ.1) THEN
IF(cros(iorba,k).EQ.1) THEN
DO lam=1, n2d
arg=ezk(lam,kz)
arg=ezk(lam,kp)
DO n=1, n2d
intw1(iorb,n)=intw1(iorb,n)+app(lam,n)*ci(iorb,lam,kz)+ &
bpp(lam,n)*di(iorb,lam,kz)
intw1(iorb,n)=intw1(iorb,n)+app(lam,n)*ci(iorb,lam,kp)+ &
bpp(lam,n)*di(iorb,lam,kp)
ENDDO
IF (ABS(DIMAG(zk(lam,kz))).LT.eps) THEN
f1p=-arg*CONJG(di(iorb,lam,kz))*zk2(lam,kz)
fl(lam,iorb)=-arg*CONJG(ci(iorb,lam,kz))*zk2(lam,kz)
IF (ABS(AIMAG(zk(lam,kp))).LT.eps) THEN
f1p=-arg*CONJG(di(iorb,lam,kp))*zk2(lam,kp)
fl(lam,iorb)=-arg*CONJG(ci(iorb,lam,kp))*zk2(lam,kp)
ELSE
f1p=-CONJG(ci(iorb,lam,kz))*zk2(lam,kz)
fl(lam,iorb)=-CONJG(di(iorb,lam,kz))*zk2(lam,kz)
f1p=-CONJG(ci(iorb,lam,kp))*zk2(lam,kp)
fl(lam,iorb)=-CONJG(di(iorb,lam,kp))*zk2(lam,kp)
ENDIF
bnlp(lam,iorb)=bnlp(lam,iorb)-f1p*emzk(lam,kz)
bnlp(lam,iorb)=bnlp(lam,iorb)-f1p*emzk(lam,kp)
ENDDO
ENDIF
ENDDO
DO iorb=1, norbnow*npol
DO iorb=1, norb*npol
iorba=iorb
IF (npol.EQ.2) iorba=(iorb+1)/2
iorbs = orbin + iorba
IF(cross(iorbs, k).EQ.1) THEN
DO iorb1=1, norbnow*npol
IF(cros(iorba,k).EQ.1) THEN
DO iorb1=1, norb*npol
iorb1a=iorb1
IF (npol.EQ.2) iorb1a=(iorb1+1)/2
iorb1s = orbin + iorb1a
tr=taunew(3,iorb1s)-rsph(nbnew(iorb1s),itnew(iorb1s))
tr=taunew(3,iorb1a)-taunew(4,iorb1a)
IF (z(k)+dz.GT.tr) THEN
c=(0.d0, 0.d0)
DO lam=1, n2d
c=c+anlp(lam,iorb1)*ci(iorb,lam,kz)+ &
bnlp(lam,iorb1)*di(iorb,lam,kz)
c=c+anlp(lam,iorb1)*ci(iorb,lam,kp)+ &
bnlp(lam,iorb1)*di(iorb,lam,kp)
ENDDO
intw2(iorb,iorb1)=intw2(iorb,iorb1)+c
ENDIF
@ -406,24 +399,25 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
!
! Rotation of linear solutions
!
CALL stop_clock('scatter')
CALL rotatef(app, bpp, bf, anlp, bnlp, bnlf, intw1, intw2, &
n2d, norbf, norbnow, npol)
!write(6,*) 'done k', k
n2d, norbf, norb, npol)
! write(6,*) 'done k', k
ENDDO
!---------------------------------------------
CALL DCOPY(2*n2d*n2d, app, 1, al, 1)
!
! To compute the 2nd half of linear solutions
!
CALL scatter_back(app, bpp, an, bn, af, ci, di, ezk, emzk, &
s1m, s2m, s3m, s4m, startk, lastk, &
norbnow, orbin, dz)
CALL scatter_back(psiper, zk, app, bpp, an, bn, af, ci, di, &
ezk, emzk, s1m, s2m, s3m, s4m, kin, kfin, &
nrz, nrzp, norb, cros, dz)
CALL DCOPY(2*n2d*n2d, bpp, 1, bl, 1)
CALL DCOPY(2*n2d*norbnow*npol, anlp, 1, anll, 1)
CALL DCOPY(2*n2d*norb*npol, anlp, 1, anll, 1)
CALL DSCAL(2*norbf*npol*2*n2d, sarea, intw1, 1)
CALL DSCAL(2*norbf*npol*norbf*npol, sarea, intw2, 1)
@ -439,58 +433,69 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
! fun1=(0.d0,0.d0)
! fund0=(0.d0,0.d0)
! fund1=(0.d0,0.d0)
k=nkofz(startk)
kz=nkofz(lastk)
DO n=1, n2d
DO lam=1, n2d
s1m(lam,n)=bf(lam,n)*ezk(lam,k)
s2m(lam,n)=al(lam,n)*ezk(lam,kz)
IF (lam.EQ.n) s2m(lam,n)=s2m(lam,n)+(1.d0,0.d0)
s3m(lam,n)=-cim*zk(lam,k)*s1m(lam,n)*tpiba
s4m(lam,n)=cim*zk(lam,kz)*s2m(lam,n)*tpiba
IF (lam.EQ.n) s4m(lam,n)=s4m(lam,n)-2.d0*cim*zk(lam,kz)*tpiba
s5m(lam,n)=bl(lam,n)*ezk(lam,k)
IF (lam.EQ.n) s5m(lam,n)=s5m(lam,n)+(1.d0,0.d0)
s6m(lam,n)=af(lam,n)*ezk(lam,kz)
s7m(lam,n)=-cim*zk(lam,k)*s5m(lam,n)*tpiba
IF (lam.EQ.n) s7m(lam,n)=s7m(lam,n)+2.d0*cim*zk(lam,k)*tpiba
s8m(lam,n)=cim*zk(lam,kz)*s6m(lam,n)*tpiba
ENDDO
ENDDO
CALL ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,k),n2d,s1m,n2d,zero,fun0,n2d)
CALL ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,kz),n2d,s2m,n2d,zero,fun1,n2d)
CALL ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,k),n2d,s3m,n2d,zero,fund0,n2d)
CALL ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,kz),n2d,s4m,n2d,zero,fund1,n2d)
CALL ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,k),n2d,s5m,n2d,zero, &
k = 1
kp = nrzp
do n=1, n2d
do lam=1, n2d
s1m(lam,n)=bf(lam,n)*ezk(lam,k)
s2m(lam,n)=al(lam,n)*ezk(lam,kp)
if (lam.eq.n) s2m(lam,n)=s2m(lam,n)+(1.d0,0.d0)
s3m(lam,n)=-cim*zk(lam,k)*s1m(lam,n)*tpiba
s4m(lam,n)=cim*zk(lam,kp)*s2m(lam,n)*tpiba
if (lam.eq.n) s4m(lam,n)=s4m(lam,n)-2.d0*cim*zk(lam,kp)*tpiba
s5m(lam,n)=bl(lam,n)*ezk(lam,k)
if (lam.eq.n) s5m(lam,n)=s5m(lam,n)+(1.d0,0.d0)
s6m(lam,n)=af(lam,n)*ezk(lam,kp)
s7m(lam,n)=-cim*zk(lam,k)*s5m(lam,n)*tpiba
if (lam.eq.n) s7m(lam,n)=s7m(lam,n)+2.d0*cim*zk(lam,k)*tpiba
s8m(lam,n)=cim*zk(lam,kp)*s6m(lam,n)*tpiba
enddo
enddo
call ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,k),n2d,s1m,n2d,zero,fun0,n2d)
call ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,kp),n2d,s2m,n2d,zero,fun1,n2d)
call ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,k),n2d,s3m,n2d,zero,fund0,n2d)
call ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,kp),n2d,s4m,n2d,zero,fund1,n2d)
call ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,k),n2d,s5m,n2d,zero, &
fun0(1,n2d+1),n2d)
CALL ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,kz),n2d,s6m,n2d,zero, &
call ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,kp),n2d,s6m,n2d,zero, &
fun1(1,n2d+1),n2d)
CALL ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,k),n2d,s7m,n2d,zero, &
call ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,k),n2d,s7m,n2d,zero, &
fund0(1,n2d+1),n2d)
CALL ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,kz),n2d,s8m,n2d,zero, &
call ZGEMM('n','n',n2d,n2d,n2d,one,psiper(1,1,kp),n2d,s8m,n2d,zero, &
fund1(1,n2d+1),n2d)
! nonlocal solutions
funl0=(0.d0,0.d0)
funl1=(0.d0,0.d0)
fundl0=(0.d0,0.d0)
fundl1=(0.d0,0.d0)
DO iorb=1, norbnow*npol
DO lam=1, n2d
do iorb = 1, norb*npol
do lam = 1, n2d
s1=ff(lam,iorb)+bnlf(lam,iorb)*ezk(lam,k)
s2=fl(lam,iorb)+anll(lam,iorb)*ezk(lam,kz)
s2=fl(lam,iorb)+anll(lam,iorb)*ezk(lam,kp)
s3=-cim*zk(lam,k)*s1*tpiba
s4=cim*zk(lam,kz)*s2*tpiba
DO ig=1, n2d
s4=cim*zk(lam,kp)*s2*tpiba
do ig=1, n2d
funl0(ig,iorb)=funl0(ig,iorb)+psiper(ig,lam,k)*s1
funl1(ig,iorb)=funl1(ig,iorb)+psiper(ig,lam,kz)*s2
funl1(ig,iorb)=funl1(ig,iorb)+psiper(ig,lam,kp)*s2
fundl0(ig,iorb)=fundl0(ig,iorb)+psiper(ig,lam,k)*s3
fundl1(ig,iorb)=fundl1(ig,iorb)+psiper(ig,lam,kz)*s4
ENDDO
ENDDO
ENDDO
fundl1(ig,iorb)=fundl1(ig,iorb)+psiper(ig,lam,kp)*s4
enddo
enddo
enddo
!
! To construct the functions in the whole rigion zin<z<zfin in the
! case of multiparallel running
!
#ifdef __PARA
CALL rotproc(fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
fundl1, intw1, intw2, n2d, norbf, norb)
#endif
CALL stop_clock('scatter_forw')
DEALLOCATE(ci)
DEALLOCATE(di)
DEALLOCATE(di)
DEALLOCATE(bf)
DEALLOCATE(an)
DEALLOCATE(bn)
@ -521,15 +526,5 @@ SUBROUTINE scatter_forw(zin, zfin, orbin, orbfin)
DEALLOCATE(emzk1)
DEALLOCATE(zk2)
!
! To construct the functions in the whole rigion zin<z<zfin in the
! case of multiparallel running
!
#ifdef __PARA
CALL rotproc(fun0, fund0, fun1, fund1, funl0, fundl0, funl1, &
fundl1, intw1, intw2, n2d, norbf, norbnow)
#endif
CALL stop_clock('scatter_forw')
RETURN
END SUBROUTINE scatter_forw

70
PWCOND/slabcpu.f90
View File

@ -1,70 +0,0 @@
!
! 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 slabcpu(nrz, nrzp, nkofz, bdl1, bdl2, bdr1, bdr2, z)
!
! To set up the division of slabs among CPU
!
USE kinds, only : DP
implicit none
integer :: k, nrz, nrzp, kin, kfin, startk, lastk
integer :: nkofz(nrz)
real(kind=DP), parameter :: eps=1.d-8
real(kind=DP) :: bdl1, bdl2, bdr1, bdr2, z(nrz)
do k=1, nrz
nkofz(k)=0
enddo
! for the left tip
do k=1, nrz
if (z(k).le.bdl1+eps) kin=k
if (z(k).le.bdl2-eps) kfin=k
enddo
if (kfin-kin.gt.0) then
call divide(kfin-kin+1,startk,lastk)
startk=startk+kin-1
lastk=lastk+kin-1
nrzp=lastk-startk+1
do k=startk, lastk
nkofz(k)=k-startk+1
enddo
endif
! for the scattering region
do k=1, nrz
if (z(k).le.bdl2+eps) kin=k
if (z(k).le.bdr1-eps) kfin=k
enddo
if (kfin-kin.gt.0) then
call divide(kfin-kin+1,startk,lastk)
startk=startk+kin-1
lastk=lastk+kin-1
do k=startk, lastk
nkofz(k)=nrzp+k-startk+1
enddo
nrzp=nrzp+lastk-startk+1
endif
! for the right tip
do k=1, nrz
if (z(k).le.bdr1+eps) kin=k
if (z(k).le.bdr2-eps) kfin=k
enddo
if (kfin-kin.gt.0) then
call divide(kfin-kin+1,startk,lastk)
startk=startk+kin-1
lastk=lastk+kin-1
do k=startk, lastk
nkofz(k)=nrzp+k-startk+1
enddo
nrzp=nrzp+lastk-startk+1
endif
return
end subroutine slabcpu

9
PWCOND/summary_tran.f90
View File

@ -5,16 +5,17 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
subroutine summary_tran(tran_file, nenergy, earr, tran_tot)
subroutine summary_tran()
!
! It writes transmission coefficients onto the file tran_file
!
USE kinds, only : DP
USE io_files, ONLY : tran_file
USE cond, ONLY : nenergy, earr, tran_tot
implicit none
character(len=80) :: tran_file
integer :: nenergy, i
real(kind=DP) :: earr(nenergy), tran_tot(nenergy)
integer :: i
!
! Output of T onto the file
!

152
PWCOND/transmit.f90
View File

@ -7,6 +7,7 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
subroutine transmit(ik, ien)
!
! This subroutine constructs the scattering states
! using the CBS of the left and right tips (saved by compbs)
@ -22,13 +23,14 @@ subroutine transmit(ik, ien)
implicit none
integer :: ik, ien, n, iorb, iorb1, iorb2, nt, &
ih, ih1, orbin, ig, ntran, ij, is, js, info
ih, ih1, ig, ntran, ij, is, js, info
integer, allocatable :: ipiv(:)
real(kind=DP) :: tk, tj, tij, eev
real(kind=DP), allocatable :: zps(:,:), eigen(:)
complex(kind=DP) :: x1, x2
complex(kind=DP), allocatable :: amat(:,:), vec1(:,:), &
tchan(:,:), veceig(:,:), zps_nc(:,:)
tmat(:,:), veceig(:,:), zps_nc(:,:), &
vec2(:,:)
eev = earr(ien)
ntran=4*n2d+npol*(norbs+nocrosl+nocrosr)
@ -40,18 +42,18 @@ implicit none
return
endif
orbin=nocrosl+noinsl
allocate( ipiv( ntran ) )
allocate( zps( norbs, norbs ) )
allocate( amat( ntran, ntran ) )
allocate( vec1( ntran, nchanl ) )
allocate( tchan( nchanr, nchanl ) )
allocate( tmat( nchanr, nchanl ) )
allocate( veceig( nchanl, nchanl ) )
allocate( eigen( nchanl ) )
if (noncolin) allocate( zps_nc( norbs*npol, norbs*npol ) )
if (noncolin) then
allocate( zps_nc( norbs*npol, norbs*npol ) )
else
allocate( zps( norbs, norbs ) )
endif
amat=(0.d0,0.d0)
!
! To form zps=zpseu-e*qq
@ -63,79 +65,26 @@ implicit none
do is=1,npol
do js=1,npol
ij=ij+1
zps_nc(npol*(iorb-1)+is, npol*(iorb1-1)+js)= &
zpseu_nc(orbin+iorb, orbin+iorb1,ij)
zps_nc(npol*(iorb-1)+is, npol*(iorb1-1)+js)= &
zpseus_nc(1,iorb,iorb1,ij)
if (lspinorb) then
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*zpseus_nc(2,iorb,iorb1,ij)
else
if (is.eq.js) &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*zpseus_nc(2,iorb,iorb1,ij)
endif
enddo
enddo
else
zps(iorb, iorb1)=zpseu(orbin+iorb,orbin+iorb1,iofspin)
zps(iorb,iorb1)=zpseus(1,iorb,iorb1)-eryd*zpseus(2,iorb,iorb1)
endif
enddo
enddo
do iorb=1, norbs-nocrosr
nt=itnew(orbin+iorb)
if(tvanp(nt)) then
do iorb1=1, norbs-nocrosr
ih=natih(orbin+iorb,2)
if (natih(orbin+iorb,1).eq.natih(orbin+iorb1,1)) then
ih1=natih(orbin+iorb1,2)
if (noncolin) then
ij=0
do is=1,npol
do js=1,npol
ij=ij+1
if (lspinorb) then
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*qq_so(ih,ih1,ij,nt)
else
if (is.eq.js) &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*qq(ih,ih1,nt)
endif
enddo
enddo
else
zps(iorb,iorb1)=zps(iorb,iorb1)- &
eryd*qq(ih,ih1,nt)
endif
endif
enddo
endif
enddo
do iorb=norbs-nocrosr+1, norbs
nt=itnew(orbin+iorb)
if(tvanp(nt)) then
ih=natih(orbin+iorb,2)
do iorb1=norbs-nocrosr+1, norbs
if (natih(orbin+iorb,1).eq.natih(orbin+iorb1,1)) then
ih1=natih(orbin+iorb1,2)
if (noncolin) then
ij=0
do is=1, npol
do js=1,npol
ij=ij+1
if (lspinorb) then
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js) &
-eryd*qq_so(ih,ih1,ij,nt)
else
if (is.eq.js) &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)= &
zps_nc(npol*(iorb-1)+is,npol*(iorb1-1)+js)- &
eryd*qq(ih,ih1,nt)
endif
enddo
enddo
else
zps(iorb,iorb1)=zps(iorb,iorb1)- &
eryd*qq(ih,ih1,nt)
endif
endif
enddo
endif
enddo
!
! Compute the part of amat which comes from the matching of
! the wave function on the boundaries
@ -251,7 +200,7 @@ implicit none
!
do n=1, nchanl
do ig=1, nchanr
tchan(ig,n)=vec1(ntran-n2d-npol*nocrosr+ig,n)
tmat(ig,n)=vec1(ntran-n2d-npol*nocrosr+ig,n)
enddo
enddo
!
@ -261,16 +210,18 @@ implicit none
do n=1, nchanl
tj=0.d0
do ig=1, nchanr
tij=DREAL(tchan(ig,n))**2+DIMAG(tchan(ig,n))**2
tij=DREAL(tmat(ig,n))**2+DIMAG(tmat(ig,n))**2
x1 = tmat(ig,n)
tj=tj+tij
WRITE( stdout,'(i5,'' --> '',i5,f12.7)') n, ig, tij
! WRITE( stdout,'(2f12.7)') dreal(x1), aimag(x1)
enddo
WRITE( stdout,'(15x,f9.5)') tj
enddo
!
! eigenchannel decomposition
!
call eigenchnl(nchanl, nchanr, tchan, veceig, eigen)
call eigenchnl(nchanl, nchanr, tmat, veceig, eigen)
WRITE( stdout,*) 'Eigenchannel decomposition:'
tk=0
do n=1, nchanl
@ -287,16 +238,53 @@ implicit none
!
tran_tot(ien) = tran_tot(ien) + wkpt(ik)*tk
!---------------------------
! Angular momentum projection of transmission
!
if(orbj_in*orbj_fin.gt.0) then
nt = orbj_fin - orbj_in + 1
allocate( vec2( ntran, nchanl ) )
x1 = (1.d0,0.d0)
x2 = (0.d0,0.d0)
call ZGEMM('n', 'n', ntran, nchanl, nchanl, x1, vec1, ntran, &
veceig, nchanl, x2, vec2, ntran)
write(6,*) 'Nchannel, Norbital, projection'
!---------------------------
! Angular momentum projection of eigenchannels
!
do n = 1, nchanl
if(eigen(n).gt.1.d-5) then
do iorb = orbj_in, orbj_fin
x1 = 0.d0
do ig = 1, 2*n2d
x1 = x1+intw1(iorb, ig)*vec2(ig, n)
enddo
do ig = 1, norbs
x1 = x1+intw2(iorb, ig)*vec2(2*n2d+ig, n)
enddo
write(6,'(2i5,f12.6)') n, iorb-orbj_in+1, DREAL(x1)**2+DIMAG(x1)**2
enddo
endif
enddo
!------------------------
deallocate(vec2)
endif
deallocate(ipiv)
deallocate(zps)
deallocate(amat)
deallocate(vec1)
deallocate(tchan)
deallocate(tmat)
deallocate(veceig)
deallocate(eigen)
if (noncolin) deallocate(zps_nc)
if (noncolin) then
deallocate(zps_nc)
else
deallocate(zps)
endif
return
end subroutine transmit