From 33c97ffce1d89540e97a5b33984079821922753f Mon Sep 17 00:00:00 2001 From: dalcorso Date: Mon, 7 Nov 2005 11:07:36 +0000 Subject: [PATCH] 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 --- PP/Makefile | 1 + PP/bands.f90 | 61 +++++++- PP/compute_sigma_avg.f90 | 304 +++++++++++++++++++++++++++++++++++++++ PW/allocate_nlpot.f90 | 1 + PW/init_us_1.f90 | 2 + 5 files changed, 363 insertions(+), 6 deletions(-) create mode 100644 PP/compute_sigma_avg.f90 diff --git a/PP/Makefile b/PP/Makefile index 93d380ac7..caff5809a 100644 --- a/PP/Makefile +++ b/PP/Makefile @@ -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 \ diff --git a/PP/bands.f90 b/PP/bands.f90 index 8d26aff6a..b3e788bf6 100644 --- a/PP/bands.f90 +++ b/PP/bands.f90 @@ -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 ! diff --git a/PP/compute_sigma_avg.f90 b/PP/compute_sigma_avg.f90 new file mode 100644 index 000000000..ebec35fac --- /dev/null +++ b/PP/compute_sigma_avg.f90 @@ -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 diff --git a/PW/allocate_nlpot.f90 b/PW/allocate_nlpot.f90 index ef10ec45c..72306a27c 100644 --- a/PW/allocate_nlpot.f90 +++ b/PW/allocate_nlpot.f90 @@ -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 diff --git a/PW/init_us_1.f90 b/PW/init_us_1.f90 index d0770c5bf..c2748f1b7 100644 --- a/PW/init_us_1.f90 +++ b/PW/init_us_1.f90 @@ -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