Commit ab6dd87f authored by Gregor Michalicek's avatar Gregor Michalicek

Initial commit for the implementation of multiple U parameters on each atom type

This is not yet tested.
parent 774881a1
...@@ -11,6 +11,7 @@ MODULE m_nmat ...@@ -11,6 +11,7 @@ MODULE m_nmat
! all atoms are stored in lda_u(), if lda_u()<0, no +U is used. ! all atoms are stored in lda_u(), if lda_u()<0, no +U is used.
! For details see Eq.(12) of Shick et al. PRB 60, 10765 (1999) ! For details see Eq.(12) of Shick et al. PRB 60, 10765 (1999)
! Part of the LDA+U package G.B., Oct. 2000 ! Part of the LDA+U package G.B., Oct. 2000
! Extension to multiple U per atom type by G.M. 2017
! ************************************************************ ! ************************************************************
CONTAINS CONTAINS
SUBROUTINE n_mat(atoms,sym, ne,usdus,jspin,we, acof,bcof,ccof, n_mmp) SUBROUTINE n_mat(atoms,sym, ne,usdus,jspin,we, acof,bcof,ccof, n_mmp)
...@@ -34,7 +35,7 @@ CONTAINS ...@@ -34,7 +35,7 @@ CONTAINS
! .. ! ..
! .. Local Scalars .. ! .. Local Scalars ..
COMPLEX c_0 COMPLEX c_0
INTEGER i,j,k,l ,mp,n,it,is,isi,natom,n_ldau,lp,m INTEGER i,j,k,l ,mp,n,it,is,isi,natom,natomTemp,n_ldau,lp,m,i_u
INTEGER ilo,ilop,ll1,nn,lmp,lm INTEGER ilo,ilop,ll1,nn,lmp,lm
REAL fac REAL fac
! .. ! ..
...@@ -45,16 +46,17 @@ CONTAINS ...@@ -45,16 +46,17 @@ CONTAINS
! !
! calculate n_mat: ! calculate n_mat:
! !
n_ldau = 0
natom = 0 natom = 0
i_u = 1
DO n = 1,atoms%ntype DO n = 1,atoms%ntype
IF (atoms%lda_u(n)%l.GE.0) THEN DO WHILE (i_u.LE.atoms%n_u)
n_ldau = n_ldau + 1 IF (atoms%lda_u(i_u)%atomType.GT.n) EXIT
n_tmp(:,:) =cmplx(0.0,0.0) natomTemp = natom
l = atoms%lda_u(n)%l n_tmp(:,:) = cmplx(0.0,0.0)
l = atoms%lda_u(i_u)%l
ll1 = (l+1)*l ll1 = (l+1)*l
DO nn = 1, atoms%neq(n) DO nn = 1, atoms%neq(n)
natom = natom + 1 natomTemp = natomTemp + 1
! !
! prepare n_mat in local frame (in noco-calculations this depends ! prepare n_mat in local frame (in noco-calculations this depends
! also on alpha(n) and beta(n) ) ! also on alpha(n) and beta(n) )
...@@ -66,8 +68,8 @@ CONTAINS ...@@ -66,8 +68,8 @@ CONTAINS
c_0 = cmplx(0.0,0.0) c_0 = cmplx(0.0,0.0)
DO i = 1,ne DO i = 1,ne
c_0 = c_0 + we(i) * ( usdus%ddn(l,n,jspin) *& c_0 = c_0 + we(i) * ( usdus%ddn(l,n,jspin) *&
conjg(bcof(i,lmp,natom))*bcof(i,lm,natom) +& conjg(bcof(i,lmp,natomTemp))*bcof(i,lm,natomTemp) +&
conjg(acof(i,lmp,natom))*acof(i,lm,natom) ) conjg(acof(i,lmp,natomTemp))*acof(i,lm,natomTemp) )
ENDDO ENDDO
n_tmp(m,mp) = c_0 n_tmp(m,mp) = c_0
ENDDO ENDDO
...@@ -85,17 +87,17 @@ CONTAINS ...@@ -85,17 +87,17 @@ CONTAINS
c_0 = cmplx(0.0,0.0) c_0 = cmplx(0.0,0.0)
DO i = 1,ne DO i = 1,ne
c_0 = c_0 + we(i) * ( usdus%uulon(ilo,n,jspin) * (& c_0 = c_0 + we(i) * ( usdus%uulon(ilo,n,jspin) * (&
conjg(acof(i,lmp,natom))*ccof(m,i,ilo,natom) +& conjg(acof(i,lmp,natomTemp))*ccof(m,i,ilo,natomTemp) +&
conjg(ccof(mp,i,ilo,natom))*acof(i,lm,natom) )& conjg(ccof(mp,i,ilo,natomTemp))*acof(i,lm,natomTemp) )&
+ usdus%dulon(ilo,n,jspin) * (& + usdus%dulon(ilo,n,jspin) * (&
conjg(bcof(i,lmp,natom))*ccof(m,i,ilo,natom) +& conjg(bcof(i,lmp,natomTemp))*ccof(m,i,ilo,natomTemp) +&
conjg(ccof(mp,i,ilo,natom))*bcof(i,lm,natom))) conjg(ccof(mp,i,ilo,natomTemp))*bcof(i,lm,natomTemp)))
ENDDO ENDDO
DO ilop = 1, atoms%nlo(n) DO ilop = 1, atoms%nlo(n)
IF (atoms%llo(ilop,n).EQ.l) THEN IF (atoms%llo(ilop,n).EQ.l) THEN
DO i = 1,ne DO i = 1,ne
c_0 = c_0 + we(i) * usdus%uloulopn(ilo,ilop,n,jspin) *& c_0 = c_0 + we(i) * usdus%uloulopn(ilo,ilop,n,jspin) *&
conjg(ccof(mp,i,ilop,natom)) *ccof(m ,i,ilo ,natom) conjg(ccof(mp,i,ilop,natomTemp)) *ccof(m ,i,ilo ,natomTemp)
ENDDO ENDDO
ENDIF ENDIF
ENDDO ENDDO
...@@ -108,10 +110,10 @@ CONTAINS ...@@ -108,10 +110,10 @@ CONTAINS
! !
! n_mmp should be rotated by D_mm' ; compare force_a21 ! n_mmp should be rotated by D_mm' ; compare force_a21
! !
DO it = 1, sym%invarind(natom) DO it = 1, sym%invarind(natomTemp)
fac = 1.0 / ( sym%invarind(natom) * atoms%neq(n) ) fac = 1.0 / ( sym%invarind(natomTemp) * atoms%neq(n) )
is = sym%invarop(natom,it) is = sym%invarop(natomTemp,it)
isi = sym%invtab(is) isi = sym%invtab(is)
d_tmp(:,:) = cmplx(0.0,0.0) d_tmp(:,:) = cmplx(0.0,0.0)
DO m = -l,l DO m = -l,l
...@@ -123,16 +125,16 @@ CONTAINS ...@@ -123,16 +125,16 @@ CONTAINS
n1_tmp = matmul( nr_tmp, d_tmp ) n1_tmp = matmul( nr_tmp, d_tmp )
DO m = -l,l DO m = -l,l
DO mp = -l,l DO mp = -l,l
n_mmp(m,mp,n_ldau) = n_mmp(m,mp,n_ldau) +conjg(n1_tmp(m,mp)) * fac n_mmp(m,mp,i_u) = n_mmp(m,mp,i_u) + conjg(n1_tmp(m,mp)) * fac
ENDDO ENDDO
ENDDO ENDDO
ENDDO ENDDO
ENDDO ! sum over equivalent atoms ENDDO ! sum over equivalent atoms
ELSE i_u = i_u + 1
natom = natom + atoms%neq(n) END DO
ENDIF natom = natom + atoms%neq(n)
ENDDO ! loop over atom types ENDDO ! loop over atom types
! do m=-l,l ! do m=-l,l
......
...@@ -229,7 +229,7 @@ CONTAINS ...@@ -229,7 +229,7 @@ CONTAINS
#endif #endif
IF (.NOT.input%secvar) THEN IF (.NOT.input%secvar) THEN
CALL timestart("hsmt extra") CALL timestart("hsmt extra")
IF (ANY(atoms%nlo>0).OR.ANY(atoms%lda_u%l.GE.0)) & IF (ANY(atoms%nlo>0).OR.(atoms%n_u.GT.0)) &
CALL hsmt_extra(DIMENSION,atoms,sym,isp,n_size,n_rank,input,nintsp,sub_comm,& CALL hsmt_extra(DIMENSION,atoms,sym,isp,n_size,n_rank,input,nintsp,sub_comm,&
hlpmsize,lmaxb,gwc,noco,l_socfirst,lapw,cell,enpara%el0,& hlpmsize,lmaxb,gwc,noco,l_socfirst,lapw,cell,enpara%el0,&
fj,gj,gk,vk,tlmplm,usdus, vs_mmp,oneD,& !in fj,gj,gk,vk,tlmplm,usdus, vs_mmp,oneD,& !in
......
...@@ -64,12 +64,11 @@ CONTAINS ...@@ -64,12 +64,11 @@ CONTAINS
COMPLEX chi11,chi21,chi22 COMPLEX chi11,chi21,chi22
INTEGER k,i,spin2,& INTEGER k,i,spin2,l,ll1,lo,jd
l,ll1,lo,jd,& INTEGER m,n,na,nn,np,i_u
m,n,na,nn,np,& INTEGER iiloh,iilos,nkvecprevath,nkvecprevats,iintsp,jintsp
iiloh,iilos,nkvecprevath,nkvecprevats,&
iintsp,jintsp
INTEGER nc,locolh,locols,nkvecprevatu,iilou,locolu INTEGER nc,locolh,locols,nkvecprevatu,iilou,locolu
INTEGER nkvecprevatuTemp,iilouTemp,locoluTemp
INTEGER ab_dim,nkvec_sv,fjstart INTEGER ab_dim,nkvec_sv,fjstart
LOGICAL enough,l_lo1 LOGICAL enough,l_lo1
! .. ! ..
...@@ -105,6 +104,7 @@ CONTAINS ...@@ -105,6 +104,7 @@ CONTAINS
na = 0 na = 0
nkvecprevats = 0 nkvecprevats = 0
nkvecprevath = 0 nkvecprevath = 0
nkvecprevatu = 0
nkvec_sv = 0 nkvec_sv = 0
!Determine index of first LO !Determine index of first LO
locols = lapw%nv(1) locols = lapw%nv(1)
...@@ -123,7 +123,10 @@ CONTAINS ...@@ -123,7 +123,10 @@ CONTAINS
iiloh = lapw%nv(1)* (lapw%nv(1)+1)/2 iiloh = lapw%nv(1)* (lapw%nv(1)+1)/2
#endif #endif
iilou = iilos
locolu = locols
i_u = 1
ntype_loop: DO n=1,atoms%ntype ntype_loop: DO n=1,atoms%ntype
IF (noco%l_noco) THEN IF (noco%l_noco) THEN
...@@ -268,29 +271,42 @@ CONTAINS ...@@ -268,29 +271,42 @@ CONTAINS
ENDIF ENDIF
END IF END IF
IF (atoms%n_u>0.and.atoms%lda_u(n)%l.GE.0.AND.gwc.EQ.1) THEN IF ((gwc.EQ.1).AND.(atoms%n_u.GT.0)) THEN
IF ( noco%l_noco .AND. (.NOT.noco%l_ss) ) THEN nkvecprevatuTemp = nkvecprevatu
CALL u_ham(& iilouTemp = iilou
atoms,input,lapw,isp,n,invsfct,& locoluTemp = locolu
ar,ai,br,bi,vs_mmp,lmaxb,& DO WHILE (i_u.LE.atoms%n_u)
alo,blo,clo,& IF (atoms%lda_u(i_u)%atomType.GT.n) EXIT
n_size,n_rank,isp,usdus,noco,& nkvecprevatuTemp = nkvecprevatu
1,1,chi11,chi22,chi21,& iilouTemp = iilou
nkvecprevatu,iilou,locolu,.false.,aa_c=aahlp) locoluTemp = locolu
ELSE IF (atoms%lda_u(i_u)%atomType.EQ.n) THEN
DO iintsp = 1,nintsp IF ((noco%l_noco).AND.(.NOT.noco%l_ss)) THEN
DO jintsp = 1,iintsp CALL u_ham(atoms,input,lapw,isp,n,i_u,invsfct,&
CALL u_ham(& ar,ai,br,bi,vs_mmp,lmaxb,&
atoms,input,lapw,isp,n,invsfct,& alo,blo,clo,&
ar,ai,br,bi,vs_mmp,lmaxb,& n_size,n_rank,isp,usdus,noco,&
alo,blo,clo,& 1,1,chi11,chi22,chi21,&
n_size,n_rank,isp,usdus,noco,& nkvecprevatuTemp,iilouTemp,locoluTemp,.false.,aa_c=aahlp)
iintsp,jintsp,chi11,chi22,chi21,& ELSE
nkvecprevatu,iilou,locolu,l_real,aa_r,aa_c) DO iintsp = 1,nintsp
ENDDO DO jintsp = 1,iintsp
ENDDO CALL u_ham(atoms,input,lapw,isp,n,i_u,invsfct,&
ENDIF ar,ai,br,bi,vs_mmp,lmaxb,&
ENDIF alo,blo,clo,&
n_size,n_rank,isp,usdus,noco,&
iintsp,jintsp,chi11,chi22,chi21,&
nkvecprevatuTemp,iilouTemp,locoluTemp,l_real,aa_r,aa_c)
END DO
END DO
END IF
END IF
i_u = i_u + 1
END DO
nkvecprevatu = nkvecprevatuTemp
iilou = iilouTemp
locolu = locoluTemp
END IF
ENDIF ! atoms%invsat(na) = 0 or 1 ENDIF ! atoms%invsat(na) = 0 or 1
!---> end loop over equivalent atoms !---> end loop over equivalent atoms
......
...@@ -14,7 +14,7 @@ MODULE m_tlmplm_store ...@@ -14,7 +14,7 @@ MODULE m_tlmplm_store
PRIVATE PRIVATE
TYPE(t_tlmplm) :: td_stored TYPE(t_tlmplm) :: td_stored
COMPLEX,ALLOCATABLE :: vs_mmp_stored(:,:,:,:) COMPLEX,ALLOCATABLE :: vs_mmp_stored(:,:,:,:)
PUBLIC write_tlmplm,read_tlmplm PUBLIC write_tlmplm, read_tlmplm, read_tlmplm_vs_mmp
CONTAINS CONTAINS
SUBROUTINE write_tlmplm(td,vs_mmp,ldau,ispin,jspin,jspins) SUBROUTINE write_tlmplm(td,vs_mmp,ldau,ispin,jspin,jspins)
TYPE(t_tlmplm),INTENT(IN) :: td TYPE(t_tlmplm),INTENT(IN) :: td
...@@ -58,15 +58,13 @@ CONTAINS ...@@ -58,15 +58,13 @@ CONTAINS
END SUBROUTINE write_tlmplm END SUBROUTINE write_tlmplm
SUBROUTINE read_tlmplm(n,jspin,nlo,ldau,tuu,tud,tdu,tdd,ind,tuulo,tuloulo,tdulo,vs_mmp) SUBROUTINE read_tlmplm(n,jspin,nlo,tuu,tud,tdu,tdd,ind,tuulo,tuloulo,tdulo)
COMPLEX,INTENT(OUT)::tuu(:),tdd(:),tud(:),tdu(:) COMPLEX,INTENT(OUT)::tuu(:),tdd(:),tud(:),tdu(:)
INTEGER,INTENT(OUT)::ind(:,:) INTEGER,INTENT(OUT)::ind(:,:)
COMPLEX,INTENT(OUT)::tuulo(:,:,:),tdulo(:,:,:),tuloulo(:,:,:) COMPLEX,INTENT(OUT)::tuulo(:,:,:),tdulo(:,:,:),tuloulo(:,:,:)
COMPLEX,INTENT(OUT)::vs_mmp(:,:)
INTEGER,INTENT(IN) :: n,jspin,nlo(:) INTEGER,INTENT(IN) :: n,jspin,nlo(:)
LOGICAL,INTENT(IN) :: ldau(:)
INTEGER:: mlo,mlolo,nn INTEGER:: mlo,mlolo
tuu=td_stored%tuu(:size(tuu,1),n,jspin) tuu=td_stored%tuu(:size(tuu,1),n,jspin)
tud=td_stored%tud(:size(tuu,1),n,jspin) tud=td_stored%tud(:size(tuu,1),n,jspin)
tdu=td_stored%tdu(:size(tuu,1),n,jspin) tdu=td_stored%tdu(:size(tuu,1),n,jspin)
...@@ -83,10 +81,18 @@ CONTAINS ...@@ -83,10 +81,18 @@ CONTAINS
tuloulo(:,:,mlolo:mlolo+nlo(n)*(nlo(n)+1)/2-1)=& tuloulo(:,:,mlolo:mlolo+nlo(n)*(nlo(n)+1)/2-1)=&
td_stored%tuloulo(:size(tuloulo,1),:size(tuloulo,2),mlolo:mlolo+nlo(n)*(nlo(n)+1)/2-1,jspin) td_stored%tuloulo(:size(tuloulo,1),:size(tuloulo,2),mlolo:mlolo+nlo(n)*(nlo(n)+1)/2-1,jspin)
ENDIF ENDIF
IF(ldau(n)) THEN
nn=count(ldau(:n-1))+1
vs_mmp=vs_mmp_stored(size(vs_mmp,1),size(vs_mmp,2),nn,jspin)
ENDIF
END SUBROUTINE read_tlmplm END SUBROUTINE read_tlmplm
SUBROUTINE read_tlmplm_vs_mmp(jspin,n_u,vs_mmp)
INTEGER, INTENT(IN) :: jspin, n_u
COMPLEX, INTENT(OUT) :: vs_mmp(:,:,:)
IF(n_u.GT.0) THEN
vs_mmp(:,:,:) = vs_mmp_stored(:,:,:,jspin)
END IF
END SUBROUTINE read_tlmplm_vs_mmp
END MODULE m_tlmplm_store END MODULE m_tlmplm_store
...@@ -57,13 +57,13 @@ CONTAINS ...@@ -57,13 +57,13 @@ CONTAINS
INTEGER, PARAMETER :: lmaxb=3 INTEGER, PARAMETER :: lmaxb=3
COMPLEX dtd,dtu,utd,utu COMPLEX dtd,dtu,utd,utu
INTEGER lo, mlotot, mlolotot, mlot_d, mlolot_d INTEGER lo, mlotot, mlolotot, mlot_d, mlolot_d
INTEGER i,ie,im,in,l1,l2,ll1,ll2,lm1,lm2,m1,m2,n,natom,m INTEGER i,ie,im,in,l1,l2,ll1,ll2,lm1,lm2,m1,m2,n,natom,m,i_u
INTEGER natrun,is,isinv,j,irinv,it INTEGER natrun,is,isinv,j,irinv,it
REAL ,PARAMETER:: zero=0.0 REAL ,PARAMETER:: zero=0.0
COMPLEX,PARAMETER:: czero=CMPLX(0.,0.) COMPLEX,PARAMETER:: czero=CMPLX(0.,0.)
! .. ! ..
! .. Local Arrays .. ! .. Local Arrays ..
COMPLEX, ALLOCATABLE :: v_mmp(:,:) COMPLEX, ALLOCATABLE :: v_mmp(:,:,:)
REAL, ALLOCATABLE :: a21(:,:),b4(:,:) REAL, ALLOCATABLE :: a21(:,:),b4(:,:)
COMPLEX forc_a21(3),forc_b4(3) COMPLEX forc_a21(3),forc_b4(3)
REAL starsum(3),starsum2(3),gvint(3),gvint2(3) REAL starsum(3),starsum2(3),gvint(3),gvint2(3)
...@@ -85,9 +85,15 @@ CONTAINS ...@@ -85,9 +85,15 @@ CONTAINS
tlmplm%tuulo(0:DIMENSION%lmd,-atoms%llod:atoms%llod,mlot_d,1),& tlmplm%tuulo(0:DIMENSION%lmd,-atoms%llod:atoms%llod,mlot_d,1),&
tlmplm%tdulo(0:DIMENSION%lmd,-atoms%llod:atoms%llod,mlot_d,1),& tlmplm%tdulo(0:DIMENSION%lmd,-atoms%llod:atoms%llod,mlot_d,1),&
tlmplm%tuloulo(-atoms%llod:atoms%llod,-atoms%llod:atoms%llod,mlolot_d,1),& tlmplm%tuloulo(-atoms%llod:atoms%llod,-atoms%llod:atoms%llod,mlolot_d,1),&
v_mmp(-lmaxb:lmaxb,-lmaxb:lmaxb),&
a21(3,atoms%nat),b4(3,atoms%nat),tlmplm%ind(0:DIMENSION%lmd,0:DIMENSION%lmd,atoms%ntype,1) ) a21(3,atoms%nat),b4(3,atoms%nat),tlmplm%ind(0:DIMENSION%lmd,0:DIMENSION%lmd,atoms%ntype,1) )
! !
IF(atoms%n_u.GT.0) THEN
ALLOCATE(v_mmp(-lmaxb:lmaxb,-lmaxb:lmaxb,atoms%n_u))
v_mmp = CMPLX(0.0,0.0)
CALL read_tlmplm_vs_mmp(jsp,atoms%n_u,v_mmp)
END IF
i_u = 1
natom = 1 natom = 1
DO n = 1,atoms%ntype DO n = 1,atoms%ntype
IF (atoms%l_geo(n)) THEN IF (atoms%l_geo(n)) THEN
...@@ -95,9 +101,9 @@ CONTAINS ...@@ -95,9 +101,9 @@ CONTAINS
forc_b4(:) = czero forc_b4(:) = czero
CALL read_tlmplm(n,jsp,atoms%nlo,atoms%lda_u%l.GE.0,& CALL read_tlmplm(n,jsp,atoms%nlo,&
tlmplm%tuu(:,n,1),tlmplm%tud(:,n,1),tlmplm%tdu(:,n,1),tlmplm%tdd(:,n,1),& tlmplm%tuu(:,n,1),tlmplm%tud(:,n,1),tlmplm%tdu(:,n,1),tlmplm%tdd(:,n,1),&
tlmplm%ind(:,:,n,1),tlmplm%tuulo(:,:,:,1),tlmplm%tuloulo(:,:,:,1),tlmplm%tdulo(:,:,:,1),v_mmp) tlmplm%ind(:,:,n,1),tlmplm%tuulo(:,:,:,1),tlmplm%tuloulo(:,:,:,1),tlmplm%tdulo(:,:,:,1))
DO natrun = natom,natom + atoms%neq(n) - 1 DO natrun = natom,natom + atoms%neq(n) - 1
a21(:,natrun) = zero a21(:,natrun) = zero
...@@ -190,9 +196,11 @@ CONTAINS ...@@ -190,9 +196,11 @@ CONTAINS
acof,bcof,ccof,aveccof,bveccof,& acof,bcof,ccof,aveccof,bveccof,&
cveccof, tlmplm,usdus, a21) cveccof, tlmplm,usdus, a21)
CALL force_a21_U(nobd,atoms,lmaxb,n,jsp,we,ne,& IF ((atoms%n_u.GT.0).AND.(i_u.LE.atoms%n_u)) THEN
usdus,v_mmp,acof,bcof,ccof,& CALL force_a21_U(nobd,atoms,lmaxb,i_u,n,jsp,we,ne,&
aveccof,bveccof,cveccof, a21) usdus,v_mmp,acof,bcof,ccof,&
aveccof,bveccof,cveccof, a21)
END IF
IF (input%l_useapw) THEN IF (input%l_useapw) THEN
! -> B4 force ! -> B4 force
DO ie = 1,ne DO ie = 1,ne
...@@ -352,7 +360,7 @@ CONTAINS ...@@ -352,7 +360,7 @@ CONTAINS
natom = natom + atoms%neq(n) natom = natom + atoms%neq(n)
ENDDO ENDDO
! !
DEALLOCATE (tlmplm%tdd,tlmplm%tuu,tlmplm%tdu,tlmplm%tud,tlmplm%tuulo,tlmplm%tdulo,tlmplm%tuloulo,v_mmp,tlmplm%ind,a21,b4) DEALLOCATE (tlmplm%tdd,tlmplm%tuu,tlmplm%tdu,tlmplm%tud,tlmplm%tuulo,tlmplm%tdulo,tlmplm%tuloulo,tlmplm%ind,a21,b4)
END SUBROUTINE force_a21 END SUBROUTINE force_a21
END MODULE m_forcea21 END MODULE m_forcea21
MODULE m_forcea21U MODULE m_forcea21U
CONTAINS CONTAINS
SUBROUTINE force_a21_U(nobd,atoms,lmaxb, itype,isp,we,ne,& SUBROUTINE force_a21_U(nobd,atoms,lmaxb,i_u,itype,isp,we,ne,&
usdus,v_mmp, acof,bcof,ccof,aveccof,bveccof,cveccof, a21) usdus,v_mmp, acof,bcof,ccof,aveccof,bveccof,cveccof, a21)
! !
!*********************************************************************** !***********************************************************************
...@@ -16,12 +16,14 @@ CONTAINS ...@@ -16,12 +16,14 @@ CONTAINS
TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_atoms),INTENT(IN) :: atoms
! .. ! ..
! .. Scalar Arguments .. ! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: nobd INTEGER, INTENT (IN) :: nobd
INTEGER, INTENT (IN) :: itype,isp,ne,lmaxb INTEGER, INTENT (IN) :: itype,isp,ne,lmaxb
INTEGER, INTENT (INOUT) :: i_u ! on input: index for the first U for atom type "itype or higher"
! on exit: index for the first U for atom type "itype+1 or higher"
! .. ! ..
! .. Array Arguments .. ! .. Array Arguments ..
REAL, INTENT (IN) :: we(nobd) REAL, INTENT (IN) :: we(nobd)
COMPLEX, INTENT (IN) :: v_mmp(-lmaxb:lmaxb,-lmaxb:lmaxb) COMPLEX, INTENT (IN) :: v_mmp(-lmaxb:lmaxb,-lmaxb:lmaxb,atoms%n_u)
COMPLEX, INTENT (IN) :: acof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat) COMPLEX, INTENT (IN) :: acof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
COMPLEX, INTENT (IN) :: bcof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat) COMPLEX, INTENT (IN) :: bcof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
COMPLEX, INTENT (IN) :: ccof(-atoms%llod:atoms%llod,nobd,atoms%nlod,atoms%nat) COMPLEX, INTENT (IN) :: ccof(-atoms%llod:atoms%llod,nobd,atoms%nlod,atoms%nat)
...@@ -42,8 +44,12 @@ CONTAINS ...@@ -42,8 +44,12 @@ CONTAINS
! comments in setlomap. ! comments in setlomap.
!*********************************************************************** !***********************************************************************
IF (atoms%lda_u(itype)%l.GE.0) THEN IF (atoms%lda_u(i_u)%atomType.GT.itype) RETURN
l = atoms%lda_u(itype)%l
DO WHILE (atoms%lda_u(i_u)%atomType.EQ.itype)
l = atoms%lda_u(i_u)%l
! !
! Add contribution for the regular LAPWs (like force_a21, but with ! Add contribution for the regular LAPWs (like force_a21, but with
! the potential matrix, v_mmp, instead of the tuu, tdd ...) ! the potential matrix, v_mmp, instead of the tuu, tdd ...)
...@@ -52,77 +58,52 @@ CONTAINS ...@@ -52,77 +58,52 @@ CONTAINS
lm = l* (l+1) + m lm = l* (l+1) + m
DO mp = -l,l DO mp = -l,l
lmp = l* (l+1) + mp lmp = l* (l+1) + mp
v_a = v_mmp(m,mp) v_a = v_mmp(m,mp,i_u)
v_b = v_mmp(m,mp) * usdus%ddn(l,itype,isp) v_b = v_mmp(m,mp,i_u) * usdus%ddn(l,itype,isp)
DO iatom = sum(atoms%neq(:itype-1))+1,sum(atoms%neq(:itype)) DO iatom = sum(atoms%neq(:itype-1))+1,sum(atoms%neq(:itype))
DO ie = 1,ne DO ie = 1,ne
DO i = 1,3 DO i = 1,3
p1 = (CONJG(acof(ie,lm,iatom)) * v_a) * aveccof(i,ie,lmp,iatom)
p2 = (CONJG(bcof(ie,lm,iatom)) * v_b) * bveccof(i,ie,lmp,iatom)
a21(i,iatom) = a21(i,iatom) + 2.0*AIMAG(p1 + p2) * we(ie)/atoms%neq(itype)
END DO
END DO
END DO
END DO ! mp
END DO ! m
p1 = ( CONJG(acof(ie,lm,iatom)) * v_a )&
* aveccof(i,ie,lmp,iatom)
p2 = ( CONJG(bcof(ie,lm,iatom)) * v_b )&
* bveccof(i,ie,lmp,iatom)
a21(i,iatom) = a21(i,iatom) + 2.0*AIMAG(&
p1 + p2 ) *we(ie)/atoms%neq(itype)
! no idea, why this did not work with ifort:
! a21(i,iatom) = a21(i,iatom) + 2.0*aimag(
! + conjg(acof(ie,lm,iatom)) * v_a *
! + *aveccof(i,ie,lmp,iatom) +
! + conjg(bcof(ie,lm,iatom)) * v_b *
! + *bveccof(i,ie,lmp,iatom) )
! + *we(ie)/neq
ENDDO
ENDDO
ENDDO
ENDDO ! mp
ENDDO ! m
! !
! If there are also LOs on this atom, with the same l as ! If there are also LOs on this atom, with the same l as
! the one of LDA+U, add another few terms ! the one of LDA+U, add another few terms
! !
DO lo = 1,atoms%nlo(itype) DO lo = 1,atoms%nlo(itype)
l = atoms%llo(lo,itype) IF (l == atoms%llo(lo,itype)) THEN
IF ( l == atoms%lda_u(itype)%l ) THEN
DO m = -l,l DO m = -l,l
lm = l* (l+1) + m lm = l* (l+1) + m
DO mp = -l,l DO mp = -l,l
lmp = l* (l+1) + mp lmp = l* (l+1) + mp
v_a = v_mmp(m,mp) v_a = v_mmp(m,mp,i_u)
v_b = v_mmp(m,mp) * usdus%uulon(lo,itype,isp) v_b = v_mmp(m,mp,i_u) * usdus%uulon(lo,itype,isp)
v_c = v_mmp(m,mp) * usdus%dulon(lo,itype,isp) v_c = v_mmp(m,mp,i_u) * usdus%dulon(lo,itype,isp)
DO iatom = sum(atoms%neq(:itype-1))+1,sum(atoms%neq(:itype)) DO iatom = sum(atoms%neq(:itype-1))+1,sum(atoms%neq(:itype))
DO ie = 1,ne DO ie = 1,ne
DO i = 1,3 DO i = 1,3
p1 = v_a * (CONJG(ccof(m,ie,lo,iatom)) * cveccof(i,mp,ie,lo,iatom))
p2 = v_b * (CONJG(acof(ie,lm,iatom)) * cveccof(i,mp,ie,lo,iatom) + &
CONJG(ccof(m,ie,lo,iatom)) * aveccof(i,ie,lmp,iatom))
p3 = v_c * (CONJG(bcof(ie,lm,iatom<