Bands.x calculates the expectation value of the spin operator on each

spinor wave-function.
(A. Smogunov and ADC)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@2427 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
dalcorso 2005-11-07 11:07:36 +00:00
parent 8796e485b6
commit 33c97ffce1
5 changed files with 363 additions and 6 deletions

View File

@ -9,6 +9,7 @@ add_shift_lc.o \
add_shift_us.o \
cft.o \
cgracsc.o \
compute_sigma_avg.o \
cube.o \
dosg.o \
do_initial_state.o \

View File

@ -20,10 +20,11 @@ PROGRAM bands
!
CHARACTER (len=256) :: filband
CHARACTER (len=256) :: outdir
LOGICAL :: lsigma(4)
INTEGER :: spin_component
INTEGER :: ios
!
NAMELIST / inputpp / outdir, prefix, filband, spin_component
NAMELIST / inputpp / outdir, prefix, filband, spin_component, lsigma
!
!
CALL start_postproc (nd_nmbr)
@ -33,6 +34,7 @@ PROGRAM bands
prefix = 'pwscf'
outdir = './'
filband = 'bands.out'
lsigma = .false.
spin_component = 1
!
IF ( npool > 1 ) CALL errore('bands','pools not implemented',npool)
@ -44,6 +46,7 @@ PROGRAM bands
READ (5, inputpp, err = 200, iostat = ios)
200 CALL errore ('do_bands', 'reading inputpp namelist', ABS (ios) )
!
lsigma(4)=.false.
tmp_dir = TRIM(outdir)
!
END IF
@ -54,6 +57,7 @@ PROGRAM bands
CALL mp_bcast( prefix, ionode_id )
CALL mp_bcast( filband, ionode_id )
CALL mp_bcast( spin_component, ionode_id )
CALL mp_bcast( lsigma(:), ionode_id )
!
! Now allocate space for pwscf variables, read and check them.
!
@ -61,14 +65,14 @@ PROGRAM bands
CALL openfil_pp
CALL init_us_1
!
CALL punch_band (filband, spin_component)
CALL punch_band (filband, spin_component, lsigma)
!
CALL stop_pp
STOP
END PROGRAM bands
!
!-----------------------------------------------------------------------
SUBROUTINE punch_band (filband, spin_component)
SUBROUTINE punch_band (filband, spin_component, lsigma)
!-----------------------------------------------------------------------
!
! This routine writes the band energies on a file. The routine orders
@ -99,6 +103,7 @@ SUBROUTINE punch_band (filband, spin_component)
COMPLEX(DP) :: pro
! the product of wavefunctions
INTEGER :: spin_component
LOGICAL :: lsigma(4)
COMPLEX(DP), ALLOCATABLE :: psiold (:,:), old (:), new (:)
! psiold: eigenfunctions at previous k-point, ordered
@ -109,7 +114,7 @@ SUBROUTINE punch_band (filband, spin_component)
COMPLEX(DP), ALLOCATABLE :: psiold_nc (:,:,:), old_nc(:,:), new_nc(:,:)
COMPLEX(DP), ALLOCATABLE :: becp_nc(:,:,:), becpold_nc(:,:,:)
! as above for the noncolinear case
INTEGER :: ibnd, jbnd, ik, ikb, ig, npwold, ios, nks1, nks2, ipol
INTEGER :: ibnd, jbnd, ik, ikb, ig, npwold, ios, nks1, nks2, ipol, ih, is1
! counters
INTEGER, ALLOCATABLE :: ok (:), igkold (:), il (:)
! ok: keeps track of which bands have been already ordered
@ -121,13 +126,25 @@ SUBROUTINE punch_band (filband, spin_component)
! ndeg : number of degenerate states
INTEGER, ALLOCATABLE :: degeneracy(:), degbands(:,:), INDEX(:)
! degbands keeps track of which states are degenerate
INTEGER :: iunpun_sigma(4)
CHARACTER(LEN=256) :: nomefile
REAL(DP), ALLOCATABLE:: edeg(:)
REAL(DP), ALLOCATABLE:: sigma_avg(:,:,:)
! expectation value of sigma
REAL(DP), PARAMETER :: eps = 0.001
! threshold (Ry) for degenerate states
COMPLEX(DP), EXTERNAL :: cgracsc, cgracsc_nc
! scalar product with the S matrix
IF (filband == ' ') RETURN
DO ipol=1,4
IF (lsigma(ipol).and..not.noncolin) THEN
CALL errore ('punch_band', 'lsigma requires noncollinear run', &
ABS (ios) )
lsigma=.false.
ENDIF
ENDDO
iunpun = 18
maxdeg = 4 * npol
!
@ -137,6 +154,17 @@ SUBROUTINE punch_band (filband, spin_component)
'formatted', err = 100, iostat = ios)
100 CALL errore ('punch_band', 'Opening filband file', ABS (ios) )
REWIND (iunpun)
DO ipol=1,4
IF (lsigma(ipol)) THEN
iunpun_sigma(ipol)=iunpun+ipol
WRITE(nomefile,'(".",i1)') ipol
OPEN (unit = iunpun_sigma(ipol), &
file = TRIM(filband)//TRIM(nomefile), &
status = 'unknown', form='formatted', err = 200, iostat = ios)
200 CALL errore ('punch_band', 'Opening filband.1 file', ABS (ios) )
REWIND (iunpun_sigma(ipol))
ENDIF
ENDDO
!
END IF
!
@ -144,6 +172,7 @@ SUBROUTINE punch_band (filband, spin_component)
ALLOCATE (psiold_nc( npwx, npol, nbnd))
ALLOCATE (becp_nc(nkb, npol, nbnd), becpold_nc(nkb, npol, nbnd))
ALLOCATE (old_nc(ngm,npol), new_nc(ngm,npol))
ALLOCATE (sigma_avg(4,nbnd,nks))
ELSE
ALLOCATE (psiold( npwx, nbnd))
ALLOCATE (old(ngm), new(ngm))
@ -193,6 +222,7 @@ SUBROUTINE punch_band (filband, spin_component)
CALL init_us_2 (npw, igk, xk (1, ik), vkb)
IF (noncolin) THEN
CALL ccalbec_nc (nkb, npwx, npw, npol, nbnd, becp_nc, vkb, evc_nc)
CALL compute_sigma_avg(sigma_avg(1,1,ik),becp_nc,ik,lsigma)
ELSE
CALL ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
END IF
@ -343,10 +373,23 @@ SUBROUTINE punch_band (filband, spin_component)
IF (ik == nks1) THEN
WRITE (iunpun, '(" &plot nbnd=",i4,", nks=",i4," /")') &
nbnd, nks2-nks1+1
DO ipol=1,4
IF (lsigma(ipol)) WRITE(iunpun_sigma(ipol), &
'(" &plot nbnd=",i4,", nks=",i4," /")') &
nbnd, nks2-nks1+1
END DO
END IF
WRITE (iunpun, '(10x,3f10.6)') xk(1,ik),xk(2,ik),xk(3,ik)
WRITE (iunpun, '(10f8.3)') (et (il (ibnd) , ik) &
WRITE (iunpun, '(10f8.3)') (et (il (ibnd) , ik) &
* rytoev, ibnd = 1, nbnd)
DO ipol=1,4
IF (lsigma(ipol)) THEN
WRITE (iunpun_sigma(ipol), '(10x,3f10.6)') &
xk(1,ik),xk(2,ik),xk(3,ik)
WRITE (iunpun_sigma(ipol), '(10f8.3)') &
(sigma_avg(ipol, il (ibnd) , ik), ibnd = 1, nbnd)
END IF
END DO
!
END IF
!
@ -357,6 +400,7 @@ SUBROUTINE punch_band (filband, spin_component)
DEALLOCATE (il, ok)
DEALLOCATE (igkold)
IF (noncolin) THEN
DEALLOCATE (sigma_avg)
DEALLOCATE (becpold_nc, becp_nc)
DEALLOCATE (new_nc, old_nc)
DEALLOCATE (psiold_nc)
@ -366,7 +410,12 @@ SUBROUTINE punch_band (filband, spin_component)
DEALLOCATE (psiold)
END IF
!
IF ( ionode ) CLOSE (iunpun)
IF ( ionode ) THEN
CLOSE (iunpun)
DO ipol=1,4
IF (lsigma(ipol)) CLOSE(iunpun_sigma(ipol))
ENDDO
ENDIF
!
RETURN
!

304
PP/compute_sigma_avg.f90 Normal file
View File

@ -0,0 +1,304 @@
!
! Copyright (C) 2005 PWSCF group
! 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 compute_sigma_avg(sigma_avg,becp_nc,ik,lsigma)
!
! This subroutine calculates the average value of the spin on
! the spinor wavefunctions.
!
USE kinds, ONLY : DP
USE noncollin_module, ONLY : noncolin, npol
USE cell_base, ONLY : alat, at, tpiba, omega
USE spin_orb, ONLY : lspinorb, fcoef
USE uspp, ONLY : nkb,qq,vkb,nhtol,nhtoj,nhtolm,indv
USE uspp_param, ONLY : nh, tvanp, nhm
USE wvfct, ONLY : nbnd, npwx, npw, igk
USE wavefunctions_module, ONLY : evc_nc, psic_nc
USE klist, ONLY : nks, xk
USE gvect, ONLY : g,gg,nr1,nr2,nr3,nrx1,nrx2,nrx3,nrxx
USE gsmooth, ONLY : nls, nlsm, nr1s, nr2s, nr3s, &
nrx1s, nrx2s, nrx3s, nrxxs, doublegrid
USE scf, ONLY : rho
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE mp_global, ONLY : me_pool
USE pffts, ONLY : npps
IMPLICIT NONE
LOGICAL :: lsigma(4)
! if true the expectation value in this direction is calculated
COMPLEX(DP) :: becp_nc(nkb,npol,nbnd)
!
REAL(KIND=DP) :: sigma_avg(4,nbnd)
INTEGER :: ik
INTEGER :: ibnd, ig, ir, ijkb0, na, np, ih, ikb, ijh, jh, jkb
INTEGER :: ipol, kh, kkb, is1, is2
INTEGER :: li, mi, lj, mj, mi1, i, j, k, ijk
REAL(DP) :: magtot1(4), magtot2(4)
REAL(DP) :: x0, y0, dx, dy, r_cut, r_aux, xx, yy
COMPLEX(DP), ALLOCATABLE :: be1(:,:), qq_lz(:,:,:)
COMPLEX(DP), ALLOCATABLE :: dfx(:), dfy(:)
COMPLEX(DP) :: c_aux, ZDOTC
IF (.NOT.(lsigma(1).OR.lsigma(2).OR.lsigma(3).OR.lsigma(4))) RETURN
ALLOCATE(be1(nhm,2))
ALLOCATE(dfx(nrxxs), dfy(nrxxs))
ALLOCATE(qq_lz(nhm,nhm,ntyp))
sigma_avg=0.d0
r_cut = 7.d0
x0 = 0.5d0*at(1,1)*alat
y0 = 0.5d0*at(2,2)*alat
dx = at(1,1)*alat/nr1s
dy = at(2,2)*alat/nr2s
qq_lz = 0.d0
DO np=1, ntyp
DO ih = 1, nh (np)
li = nhtol(ih,np)
mi = nhtolm(ih,np) - li**2
IF (mi.EQ.2) THEN
mi1 = 3
c_aux = -(0.d0,1.d0)
ELSE IF (mi.EQ.3) THEN
mi1 = 2
c_aux = (0.d0,1.d0)
ELSE IF (mi.EQ.4) THEN
mi1 = 5
c_aux = -(0.d0,2.d0)
ELSE IF (mi.EQ.5) THEN
mi1 = 4
c_aux = (0.d0,2.d0)
END IF
DO jh = ih+1, nh (np)
lj = nhtol(jh,np)
mj = nhtolm(jh,np) - lj**2
IF (lj.EQ.li.AND.mj.EQ.mi1) THEN
IF (mj.GT.mi) THEN
r_aux = qq(ih,jh-1,np)
ELSE
r_aux = qq(ih,jh+1,np)
END IF
qq_lz(ih,jh,np) = c_aux * r_aux
END IF
END DO
END DO
DO ih = 1, nh (np)
DO jh = 1, ih-1
qq_lz(ih,jh,np) = CONJG(qq_lz(jh,ih,np))
END DO
END DO
END DO
DO ibnd = 1, nbnd
rho = 0.d0
magtot1 = 0.d0
magtot2 = 0.d0
!-- Pseudo part
psic_nc = (0.D0,0.D0)
DO ipol=1,npol
DO ig = 1, npw
psic_nc(nls(igk(ig)),ipol)=evc_nc(ig,ipol,ibnd)
END DO
call cft3s (psic_nc(1,ipol), nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
ENDDO
!
! Calculate the three components of the magnetization
! (stored in rho(ir,2-4) )
!
IF (lsigma(1)) THEN
DO ir = 1,nrxxs
rho(ir,2) = rho(ir,2) + 2.D0* &
(real(psic_nc(ir,1))*real(psic_nc(ir,2)) + &
DIMAG(psic_nc(ir,1))*DIMAG(psic_nc(ir,2)))
END DO
IF (doublegrid) CALL interpolate( rho(1,2), rho(1,2), 1 )
END IF
IF (lsigma(2)) THEN
DO ir = 1,nrxxs
rho(ir,3) = rho(ir,3) + 2.D0* &
(real(psic_nc(ir,1))*DIMAG(psic_nc(ir,2)) - &
real(psic_nc(ir,2))*DIMAG(psic_nc(ir,1)))
END DO
IF (doublegrid) CALL interpolate( rho(1,3), rho(1,3), 1 )
END IF
IF (lsigma(3)) THEN
DO ir = 1,nrxxs
rho(ir,4) = rho(ir,4) + &
(real(psic_nc(ir,1))**2+DIMAG(psic_nc(ir,1))**2 &
-real(psic_nc(ir,2))**2-DIMAG(psic_nc(ir,2))**2)
END DO
IF (doublegrid) CALL interpolate( rho(1,4), rho(1,4), 1 )
END IF
IF (lsigma(4)) THEN
!-- Calculate pseudo part of L_z
DO ipol = 1, npol
dfx = 0.d0
dfy = 0.d0
dfx(nls(igk(1:npw))) = (xk(1,ik)+g(1,igk(1:npw)))*tpiba* &
(0.d0,1.d0)*evc_nc(1:npw,ipol,ibnd)
dfy(nls(igk(1:npw))) = (xk(2,ik)+g(2,igk(1:npw)))*tpiba* &
(0.d0,1.d0)*evc_nc(1:npw,ipol,ibnd)
CALL cft3s( dfx, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
CALL cft3s( dfy, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2 )
DO i = 1, nr1s
xx = (i-1)*dx - x0
DO j = 1, nr2s
yy = (j-1)*dy - y0
r_aux = SQRT (xx**2 + yy**2)
IF (r_aux.LE.r_cut) THEN
DO k = 1, npps(me_pool+1)
ijk = i + (j-1) * nrx1s + (k-1) * nrx1s * nrx2s
dfx(ijk) = xx * dfy(ijk) - yy * dfx(ijk)
END DO
ELSE
DO k = 1, npps(me_pool+1)
ijk = i + (j-1) * nrx1s + (k-1) * nrx1s * nrx2s
dfx (ijk) = 0.d0
END DO
END IF
END DO
END DO
c_aux = ZDOTC(nrxxs, psic_nc(1,ipol), 1, dfx, 1)
magtot1(4) = magtot1(4) + DIMAG(c_aux)
END DO
CALL reduce( 1, magtot1(4) )
magtot1(4) = magtot1(4)/(nr1s*nr2s*nr3s)
END IF
DO ipol=1,3
IF (lsigma(ipol)) THEN
DO ir = 1,nrxx
magtot1(ipol) = magtot1(ipol) + rho(ir,ipol+1)
END DO
CALL reduce( 1, magtot1(ipol) )
magtot1(ipol) = magtot1(ipol) / ( nr1 * nr2 * nr3 )
END IF
END DO
!-- Augmentation part
ijkb0 = 0
!
DO np = 1, ntyp
!
IF ( tvanp(np) ) THEN
!
DO na = 1, nat
!
IF (ityp(na)==np) THEN
!
be1 = 0.d0
DO ih = 1, nh(np)
ikb = ijkb0 + ih
IF (lspinorb) THEN
DO kh = 1, nh(np)
IF ((nhtol(kh,np)==nhtol(ih,np)).AND. &
(nhtoj(kh,np)==nhtoj(ih,np)).AND. &
(indv(kh,np)==indv(ih,np))) THEN
kkb=ijkb0 + kh
DO is1=1,2
DO is2=1,2
be1(ih,is1)=be1(ih,is1)+ &
fcoef(ih,kh,is1,is2,np)* &
becp_nc(kkb,is2,ibnd)
END DO
END DO
END IF
END DO
ELSE
DO is1=1,2
be1(ih,is1) = becp_nc(ikb,is1,ibnd)
END DO
END IF
END DO
IF (lsigma(1)) THEN
DO ih = 1, nh(np)
magtot2(1)=magtot2(1)+ 2.d0*qq(ih,ih,np)*DREAL &
( be1(ih,2)*DCONJG(be1(ih,1)) )
DO jh = ih + 1, nh(np)
magtot2(1)=magtot2(1)+2.d0*qq(ih,jh,np)*DREAL &
(be1(jh,2)*DCONJG(be1(ih,1))+ &
be1(jh,1)*DCONJG(be1(ih,2)) )
ENDDO
ENDDO
ENDIF
IF (lsigma(2)) THEN
DO ih = 1, nh(np)
magtot2(2)=magtot2(2)+ 2.d0*qq(ih,ih,np)*DIMAG &
( be1(ih,2)*DCONJG(be1(ih,1)) )
DO jh = ih + 1, nh(np)
magtot2(2)=magtot2(2) + 2.d0*qq(ih,jh,np)*DIMAG &
( be1(jh,2) * DCONJG(be1(ih,1)) &
- be1(jh,1) * DCONJG(be1(ih,2)) )
END DO
END DO
END IF
IF (lsigma(3)) THEN
DO ih = 1, nh(np)
magtot2(3) = magtot2(3) + qq(ih,ih,np)* &
( ABS(be1(ih,1))**2 - ABS(be1(ih,2))**2 )
DO jh = ih + 1, nh(np)
magtot2(3) = magtot2(3) + 2.d0*qq(ih,jh,np)*DREAL &
(be1(jh,1)*DCONJG(be1(ih,1)) &
-be1(jh,2)*DCONJG(be1(ih,2)) )
END DO
END DO
END IF
IF (lsigma(4)) THEN
DO ih = 1, nh(np)
DO jh = ih + 1, nh(np)
magtot2(4)= magtot2(4)+2.d0*DREAL(qq_lz(ih,jh,np)* &
( CONJG(be1(ih,1))*be1(jh,1) + &
CONJG(be1(ih,2))*be1(jh,2) ) )
END DO
END DO
END IF
!
ijkb0 = ijkb0 + nh(np)
!
END IF
!
END DO
!
ELSE
!
DO na = 1, nat
!
IF ( ityp(na) == np ) ijkb0 = ijkb0 + nh(np)
!
END DO
!
END IF
!
END DO
DO ipol=1,3
IF (lsigma(ipol)) &
sigma_avg(ipol,ibnd) = 0.5d0 * ( magtot1(ipol) + magtot2(ipol) )
END DO
IF (lsigma(4)) &
sigma_avg(4,ibnd) = magtot1(4) + magtot2(4) + sigma_avg(3,ibnd)
END DO
DEALLOCATE(be1)
DEALLOCATE(dfx,dfy)
DEALLOCATE(qq_lz)
RETURN
END SUBROUTINE compute_sigma_avg

View File

@ -95,6 +95,7 @@ subroutine allocate_nlpot
endif
if (lspinorb) then
allocate (qq_so(nhm, nhm, 4, nsp))
allocate (qq( nhm, nhm, nsp))
allocate (dvan_so( nhm, nhm, nspin, nsp))
allocate (fcoef(nhm,nhm,2,2,nsp))
else

View File

@ -287,6 +287,7 @@ subroutine init_us_1
do ih=1,nh(nt)
do jh=1,nh(nt)
call qvan2 (1, ih, jh, nt, gg, qgm, ylmk0)
qq (ih, jh, nt) = omega * DBLE (qgm (1) )
do kh=1,nh(nt)
do lh=1,nh(nt)
ijs=0
@ -319,6 +320,7 @@ subroutine init_us_1
100 continue
if (lspinorb) then
call reduce ( nhm * nhm * ntyp * 8, qq_so )
call reduce ( nhm * nhm * ntyp, qq )
else
call reduce ( nhm * nhm * ntyp, qq )
endif