diff --git a/D3/d0rhod2v.f90 b/D3/d0rhod2v.f90 index 839fc9bb8..61b21e30a 100644 --- a/D3/d0rhod2v.f90 +++ b/D3/d0rhod2v.f90 @@ -173,7 +173,7 @@ subroutine d0rhod2v (ipert, drhoscf) d3dywrk (na_icart, na_jcart) = d3dywrk (na_icart, na_jcart) & + (alpha(1)*alpha(2) + alpha(3)*alpha(4) & - alpha(5)*alpha(6) - alpha(7)*alpha(8)) * & - dvan (ikb,ikb,nt) * wgg * 2.0d0 + dvan (ikb,ikb,1,nt) * wgg * 2.0d0 enddo end if enddo diff --git a/D3/d3vrho.f90 b/D3/d3vrho.f90 index cc6c4a6ba..4fa124990 100644 --- a/D3/d3vrho.f90 +++ b/D3/d3vrho.f90 @@ -129,7 +129,7 @@ subroutine d3vrho call reduce (16, alpha) #endif d3dynwrk (na_k, na_i, na_j) = d3dynwrk (na_k, na_i, na_j) - & - 2.0d0 * dvan(ikb,ikb,nt) * wgg * & + 2.0d0 * dvan(ikb,ikb,1,nt) * wgg * & DIMAG(alpha(1)*alpha(2) + alpha(3)*alpha(4) + & alpha(5)*alpha(6) + alpha(7)*alpha(8)) enddo diff --git a/D3/dqrhod2v.f90 b/D3/dqrhod2v.f90 index 00f61d814..54ffca19f 100644 --- a/D3/dqrhod2v.f90 +++ b/D3/dqrhod2v.f90 @@ -167,7 +167,7 @@ subroutine dqrhod2v (ipert, drhoscf) d3dywrk(na_icart,na_jcart) = d3dywrk(na_icart,na_jcart) & + conjg(alpha(1) * alpha(2) + alpha(3) * alpha(4) - & alpha(5) * alpha(6) - alpha(7) * alpha(8) ) & - * dvan (ikb, ikb, nt) * wgg * 2.0d0 + * dvan (ikb, ikb, 1, nt) * wgg * 2.0d0 enddo endif enddo diff --git a/D3/dvdpsi.f90 b/D3/dvdpsi.f90 index 21fe90b0f..c1a55050e 100644 --- a/D3/dvdpsi.f90 +++ b/D3/dvdpsi.f90 @@ -88,9 +88,9 @@ subroutine dvdpsi (nu_i, xq_, dvloc, vkb_, vkbq_, psi_, dvpsi_) g (3, igk (ig) ) * u_x (mu + 3, nu_i) ) ) enddo do ibnd = 1, nbnd - ps(1,ibnd) = dvan(ikb,ikb,nt) * & + ps(1,ibnd) = dvan(ikb,ikb,1,nt) * & ZDOTC(npw, wrk2, 1, psi_(1,ibnd), 1) - ps(2,ibnd) = dvan(ikb,ikb,nt) * & + ps(2,ibnd) = dvan(ikb,ikb,1,nt) * & ZDOTC(npw,vkb_(1,jkb),1,psi_(1,ibnd),1) enddo #ifdef __PARA diff --git a/Gamma/dvpsi_e.f90 b/Gamma/dvpsi_e.f90 index 3a7be852f..fe15a1abe 100644 --- a/Gamma/dvpsi_e.f90 +++ b/Gamma/dvpsi_e.f90 @@ -87,8 +87,8 @@ subroutine dvpsi_e(kpoint,ipol) do ih=1,nh(nt) jkb=jkb+1 do ibnd = 1,nbnd - dbec_(jkb,ibnd) = dbec(jkb,ibnd)*dvan(ih,ih,nt) - becp_(jkb,ibnd) = becp(jkb,ibnd)*dvan(ih,ih,nt) + dbec_(jkb,ibnd) = dbec(jkb,ibnd)*dvan(ih,ih,1,nt) + becp_(jkb,ibnd) = becp(jkb,ibnd)*dvan(ih,ih,1,nt) enddo end do end if diff --git a/Gamma/dvpsi_kb.f90 b/Gamma/dvpsi_kb.f90 index ae1b538db..ed57295ce 100644 --- a/Gamma/dvpsi_kb.f90 +++ b/Gamma/dvpsi_kb.f90 @@ -114,8 +114,8 @@ subroutine dvpsi_kb(kpoint,nu) ! do ibnd = 1,nbnd do ih = 1,nh(nt) - bec1(ih,ibnd) = dvan(ih,ih,nt) * bec1(ih,ibnd) - bec2(ih,ibnd) = dvan(ih,ih,nt) * bec2(ih,ibnd) + bec1(ih,ibnd) = dvan(ih,ih,1,nt) * bec1(ih,ibnd) + bec2(ih,ibnd) = dvan(ih,ih,1,nt) * bec2(ih,ibnd) end do end do ! diff --git a/Gamma/rhod2vkb.f90 b/Gamma/rhod2vkb.f90 index 363e11fd0..02fe9eeb6 100644 --- a/Gamma/rhod2vkb.f90 +++ b/Gamma/rhod2vkb.f90 @@ -147,7 +147,7 @@ subroutine rhod2vkb(dyn0) dynkb(nu_i,nu_j) = dynkb(nu_i,nu_j) + & (-becp1(jkb+ih,ibnd,ipol)*becp1(jkb+ih,ibnd,jpol) & +becp2(jkb+ih,ibnd,ijpol)*becp(jkb+ih,ibnd) ) & - * dvan(ih,ih,nt) * weight + * dvan(ih,ih,1,nt) * weight end do end do end do diff --git a/Modules/atom.f90 b/Modules/atom.f90 index 3273e600d..2950ab380 100644 --- a/Modules/atom.f90 +++ b/Modules/atom.f90 @@ -22,6 +22,7 @@ MODULE atom dx(npsx), &! linear interval for logaritmic mesh r(ndmx,npsx), &! radial logaritmic mesh rab(ndmx,npsx), &! derivative of the radial mesh + jchi(nchix,npsx), &! total angular momentum of atomic orbitals chi(ndmx,nchix,npsx), &! radial atomic orbitals oc(nchix,npsx), &! atomic level occupation rho_at(ndmx,npsx), &! radial atomic charge density diff --git a/PP/pw2casino.f90 b/PP/pw2casino.f90 index 1bb163337..480db1abf 100644 --- a/PP/pw2casino.f90 +++ b/PP/pw2casino.f90 @@ -172,13 +172,13 @@ subroutine compute_casino do ih = 1, nh (nt) ikb = ijkb0 + ih enl=enl+conjg(becp(ikb,ibnd))*becp(ikb,ibnd) & - *wg(ibnd,ikk)* dvan(ih,ih,nt) + *wg(ibnd,ikk)* dvan(ih,ih,1,nt) DO jh = ( ih + 1 ), nh(nt) jkb = ijkb0 + jh enl=enl + & (conjg(becp(ikb,ibnd))*becp(jkb,ibnd)+& conjg(becp(jkb,ibnd))*becp(ikb,ibnd))& - * wg(ibnd,ikk) * dvan(ih,jh,nt) + * wg(ibnd,ikk) * dvan(ih,jh,1,nt) END DO diff --git a/PW/addusforce.f90 b/PW/addusforce.f90 index 7ddb01e76..4b1895dd6 100644 --- a/PW/addusforce.f90 +++ b/PW/addusforce.f90 @@ -119,7 +119,7 @@ subroutine addusforce (forcenl) ! end do ! WRITE( stdout,'( "dion pseudo ",i4)') nt ! do ih = 1, nh(nt) - ! WRITE( stdout,'(8f9.4)') (dvan(ih,jh,nt),jh=1,nh(nt)) + ! WRITE( stdout,'(8f9.4)') (dvan(ih,jh,1,nt),jh=1,nh(nt)) ! end do do is = 1, nspin do na = 1, nat diff --git a/PW/allocate_nlpot.f90 b/PW/allocate_nlpot.f90 index 6b17a6ff0..667ba1d7f 100644 --- a/PW/allocate_nlpot.f90 +++ b/PW/allocate_nlpot.f90 @@ -35,7 +35,8 @@ subroutine allocate_nlpot USE wvfct, ONLY: npwx, npw, igk, igk_l2g, g2kin USE us, ONLY: nh, indv, nhtol, nhtolm, qq, dvan, deeq, qrad, vkb, tab, & tab_at, dq, qgm, becsum, nhm, lqx, nqx, nqxq, nkb, lmaxkb, lll, nbeta, & - tvanp, newpseudo + tvanp, newpseudo, nhtoj + USE spin_orb, ONLY: lspinorb, qq_spinorb, fcoef implicit none ! ! a few local variables @@ -90,12 +91,18 @@ subroutine allocate_nlpot nkb = nkb + nh (ityp(na)) enddo ! - allocate (indv( nhm, ntyp)) - allocate (nhtol( nhm, ntyp)) + allocate (indv( nhm, ntyp)) + allocate (nhtol(nhm, ntyp)) allocate (nhtolm(nhm, ntyp)) + allocate (nhtoj(nhm, ntyp)) allocate (qq( nhm, nhm, ntyp)) - allocate (dvan( nhm, nhm, ntyp)) + allocate (dvan( nhm, nhm, nspin, ntyp)) allocate (deeq( nhm, nhm, nat, nspin)) + if (lspinorb) then + allocate (qq_spinorb(nhm, nhm, 4, ntyp)) + allocate (fcoef(nhm,nhm,2,2,ntyp)) + endif + ! nqxq = ( (sqrt(gcutm) + sqrt(xqq(1)**2 + xqq(2)**2 + xqq(3)**2) ) & / dq + 4) * cell_factor diff --git a/PW/clean_pw.f90 b/PW/clean_pw.f90 index df6acfcc0..a07accb66 100644 --- a/PW/clean_pw.f90 +++ b/PW/clean_pw.f90 @@ -25,7 +25,7 @@ SUBROUTINE clean_pw() USE relax, ONLY : if_pos USE wavefunctions_module, ONLY : evc, psic USE us, ONLY : indv, nhtol, nhtolm, qq, dvan, deeq, qrad, & - vkb, qgm, becsum, tab, tab_at + vkb, qgm, becsum, tab, tab_at, nhtoj USE ldaU, ONLY : ns, nsnew, swfcatom USE extfield, ONLY : forcefield USE becmod, ONLY : becp @@ -36,7 +36,9 @@ SUBROUTINE clean_pw() USE afftnec, ONLY : auxp #endif USE fft_types, ONLY : fft_dlay_deallocate + USE spin_orb, ONLY : lspinorb, qq_spinorb, fcoef USE constraints_module, ONLY : constr, target + ! IMPLICIT NONE ! @@ -96,7 +98,8 @@ SUBROUTINE clean_pw() IF ( ALLOCATED( g2kin ) ) DEALLOCATE( g2kin ) IF ( ALLOCATED( indv ) ) DEALLOCATE( indv ) IF ( ALLOCATED( nhtol ) ) DEALLOCATE( nhtol ) - IF ( ALLOCATED( nhtolm) ) DEALLOCATE( nhtolm) + IF ( ALLOCATED( nhtolm ) ) DEALLOCATE( nhtolm ) + IF ( ALLOCATED( nhtoj ) ) DEALLOCATE( nhtoj ) IF ( ALLOCATED( qq ) ) DEALLOCATE( qq ) IF ( ALLOCATED( dvan ) ) DEALLOCATE( dvan ) IF ( ALLOCATED( deeq ) ) DEALLOCATE( deeq ) @@ -108,6 +111,11 @@ SUBROUTINE clean_pw() IF ( ALLOCATED( nsnew ) ) DEALLOCATE( nsnew ) IF ( ALLOCATED( tab ) ) DEALLOCATE( tab ) IF ( ALLOCATED( tab_at ) ) DEALLOCATE( tab_at ) + IF (lspinorb) then + IF ( ALLOCATED( qq_spinorb ) ) DEALLOCATE( qq_spinorb ) + IF ( ALLOCATED( fcoef ) ) DEALLOCATE( fcoef ) + END IF + ! ! ... arrays allocated in allocate_wfc.f90 ( and never deallocated ) ! diff --git a/PW/init_us_1.f90 b/PW/init_us_1.f90 index 086564ba9..fbb14c072 100644 --- a/PW/init_us_1.f90 +++ b/PW/init_us_1.f90 @@ -70,7 +70,7 @@ subroutine init_us_1 allocate (besr( ndm)) allocate (qtot( ndm , nbrx , nbrx)) allocate (ylmk0( lqx * lqx)) - dvan (:,:,:) = 0.d0 + dvan (:,:,:,:) = 0.d0 qq (:,:,:) = 0.d0 ap (:,:,:) = 0.d0 if (lqx > 0) qrad(:,:,:,:)= 0.d0 @@ -104,7 +104,7 @@ subroutine init_us_1 nhtolm(ih, nt) == nhtolm(jh, nt) ) then ir = indv (ih, nt) is = indv (jh, nt) - dvan (ih, jh, nt) = dion (ir, is, nt) + dvan (ih, jh, 1, nt) = dion (ir, is, nt) endif enddo enddo diff --git a/PW/newd.f90 b/PW/newd.f90 index a69130001..44df7a7db 100644 --- a/PW/newd.f90 +++ b/PW/newd.f90 @@ -42,7 +42,7 @@ subroutine newd do is = 1, nspin do ih = 1, nh (nt) do jh = ih, nh (nt) - deeq (ih, jh, na, is) = dvan (ih, jh, nt) + deeq (ih, jh, na, is) = dvan (ih, jh, 1, nt) deeq (jh, ih, na, is) = deeq (ih, jh, na, is) enddo enddo @@ -149,14 +149,14 @@ subroutine newd ! end do do ih = 1, nh (nt) do jh = ih, nh (nt) - deeq (ih, jh, na, is) = deeq (ih, jh, na, is) + dvan (ih, jh, nt) + deeq (ih, jh, na, is) = deeq (ih, jh, na, is) + dvan (ih,jh,1,nt) deeq (jh, ih, na, is) = deeq (ih, jh, na, is) enddo enddo enddo ! WRITE( stdout,'( "dion pseudo ",i4)') nt ! do ih = 1, nh(nt) - ! WRITE( stdout,'(8f9.4)') (dvan(ih,jh,nt),jh=1,nh(nt)) + ! WRITE( stdout,'(8f9.4)') (dvan(ih,jh,1,nt),jh=1,nh(nt)) ! end do enddo diff --git a/PW/pwcom.f90 b/PW/pwcom.f90 index 8080740b3..5c72dcb15 100644 --- a/PW/pwcom.f90 +++ b/PW/pwcom.f90 @@ -441,6 +441,8 @@ MODULE us dion(nbrx,nbrx,npsx), &! D_{mu,nu} parameters (in the ! atomic case) betar(ndmx,nbrx,npsx), &! radial beta_{mu} functions + jjj(nbrx,npsx), &! total angular momentum of + ! the beta function qqq(nbrx,nbrx,npsx), &! q_{mu,nu} parameters (in the ! atomic case) qfunc(ndmx,nbrx,nbrx,npsx), &! Q_{mu,nu}(|r|) function for @@ -471,14 +473,15 @@ MODULE us nhtolm(:,:) ! correspondence n <-> combined lm index for (l,m) COMPLEX(KIND=DP), ALLOCATABLE, TARGET :: & vkb(:,:), &! all beta functions in reciprocal space + dvan(:,:,:,:), &! the D functions of the solid + deeq(:,:,:,:), &! the integral of V_eff and Q_{nm} qgm(:) ! complete fourier transform of Q REAL(KIND=DP), PARAMETER:: & dq = 0.01D0 ! space between points in the pseudopotential tab. REAL(KIND=DP), ALLOCATABLE :: & qq(:,:,:), &! the q functions in the solid - dvan(:,:,:), &! the D functions of the solid - deeq(:,:,:,:), &! the integral of V_eff and Q_{nm} becsum(:,:,:), &! the sum of bec functions + nhtoj(:,:), &! correspondence n <-> total angular momentum qrad(:,:,:,:), &! radial FT of Q functions tab(:,:,:), &! interpolation table for PPs tab_at(:,:,:) ! interpolation table for atomic wfc @@ -596,6 +599,29 @@ MODULE fixed_occ tfixed_occ ! if .TRUE. the occupations are fixed. ! END MODULE fixed_occ + +MODULE spin_orb + + USE kinds, ONLY: DP + USE parameters, ONLY : lmaxx + + SAVE + + LOGICAL :: & + lspinorb ! if .TRUE. this is a spin-robit calculation + + COMPLEX (kind=dp) :: rot_ylm(2*lmaxx+1,2*lmaxx+1) ! transform real + ! spherical harmonics into complex ones + COMPLEX (kind=dp), ALLOCATABLE :: fcoef(:,:,:,:,:) ! function needed to + ! account for spinors. + COMPLEX (kind=dp),ALLOCATABLE :: qq_spinorb(:,:,:,:) ! the qq coefficients + ! in the spin_orbit case. NB: This variable should + ! not be necessary, but to introduce an additional + ! index in any place where qq is used was too + ! cumbersome. As a first step I have introduced + ! another variable with the required index + ! (ADC 26/3/04) +END MODULE spin_orb ! ! MODULE pwcom @@ -626,5 +652,6 @@ MODULE pwcom USE sticks USE bp USE fixed_occ + USE spin_orb ! END MODULE pwcom