Commit 9d9906c2 authored by Daniel Wortmann's avatar Daniel Wortmann

Removed jij data-type, added first implementation of further force-theorem...

Removed jij data-type, added first implementation of further force-theorem modes, these lack tests&support in xml-shema
parent 6ddaa944
......@@ -8,7 +8,7 @@ MODULE m_genNewNocoInp
CONTAINS
SUBROUTINE genNewNocoInp(input,atoms,jij,noco,noco_new)
SUBROUTINE genNewNocoInp(input,atoms,noco,noco_new)
USE m_juDFT
USE m_types
......@@ -19,7 +19,6 @@ SUBROUTINE genNewNocoInp(input,atoms,jij,noco,noco_new)
TYPE(t_input),INTENT(IN) :: input
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_jij),INTENT(IN) :: jij
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_noco),INTENT(INOUT) :: noco_new
......@@ -50,7 +49,7 @@ SUBROUTINE genNewNocoInp(input,atoms,jij,noco,noco_new)
OPEN (24,file='nocoinp',form='formatted', status='old')
REWIND (24)
CALL rw_noco_write(atoms,jij,noco_new, input)
CALL rw_noco_write(atoms,noco_new, input)
CLOSE (24)
END SUBROUTINE genNewNocoInp
......
......@@ -19,7 +19,7 @@ CONTAINS
!> The matrices generated and diagonalized here are of type m_mat as defined in m_types_mat.
!>@author D. Wortmann
SUBROUTINE eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
sym,kpts,DIMENSION, vacuum, input, cell, enpara_in,enpara,banddos, noco,jij, oneD,hybrid,&
sym,kpts,DIMENSION, vacuum, input, cell, enpara_in,enpara,banddos, noco, oneD,hybrid,&
it,eig_id,results,inden,v,vx)
USE m_constants, ONLY : pi_const,sfp_const
USE m_types
......@@ -55,7 +55,6 @@ CONTAINS
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_banddos),INTENT(IN) :: banddos
TYPE(t_jij),INTENT(IN) :: jij
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
......@@ -154,7 +153,7 @@ CONTAINS
CALL lapw%init(input,noco, kpts,atoms,sym,nk,cell,l_zref, mpi)
call timestart("Setup of H&S matrices")
CALL eigen_hssetup(jsp,mpi,DIMENSION,hybrid,enpara,input,vacuum,noco,jij,sym,&
CALL eigen_hssetup(jsp,mpi,DIMENSION,hybrid,enpara,input,vacuum,noco,sym,&
stars,cell,sphhar,atoms,ud,td,v,lapw,l_real,smat,hmat)
CALL timestop("Setup of H&S matrices")
......@@ -191,13 +190,8 @@ CONTAINS
#else
ne_found=ne_all
#endif
!jij%eig_l = 0.0 ! need not be used, if hdf-file is present
IF (.NOT.l_real) THEN
IF (.NOT.jij%l_J) THEN
zMat%data_c(:lapw%nmat,:ne_found) = CONJG(zMat%data_c(:lapw%nmat,:ne_found))
ELSE
zMat%data_c(:lapw%nmat,:ne_found) = CMPLX(0.0,0.0)
ENDIF
ENDIF
CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,&
eig(:ne_found),n_start=mpi%n_size,n_end=mpi%n_rank,zmat=zMat)
......
......@@ -15,7 +15,7 @@ CONTAINS
!! 4. The vacuum part is added (in hsvac())
!! 5. The matrices are copied to the final matrix, in the noco-case the full matrix is constructed from the 4-parts.
SUBROUTINE eigen_hssetup(isp,mpi,DIMENSION,hybrid,enpara,input,vacuum,noco,jij,sym,&
SUBROUTINE eigen_hssetup(isp,mpi,DIMENSION,hybrid,enpara,input,vacuum,noco,sym,&
stars,cell,sphhar,atoms,ud,td,v,lapw,l_real,smat_final,hmat_final)
USE m_hs_int
USE m_hsvac
......@@ -33,7 +33,6 @@ CONTAINS
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_jij),INTENT(IN) :: jij
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
......@@ -82,7 +81,7 @@ CONTAINS
IF (input%film) THEN
CALL timestart("Vacuum part")
CALL hsvac(vacuum,stars,DIMENSION, atoms,mpi,isp,input,v,enpara%evac0,cell,&
lapw,sym, noco,jij,hmat,smat)
lapw,sym, noco,hmat,smat)
CALL timestop("Vacuum part")
ENDIF
!Now copy the data into final matrix
......
......@@ -12,7 +12,7 @@ CONTAINS
!-----------------------------------------------------------
SUBROUTINE hsvac(&
vacuum,stars,DIMENSION, atoms,mpi,jsp,input,v,evac,cell,&
lapw,sym, noco,jij,hmat,smat)
lapw,sym, noco,hmat,smat)
USE m_vacfun
......@@ -22,7 +22,6 @@ CONTAINS
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_jij),INTENT(IN) :: jij
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
......@@ -57,13 +56,10 @@ CONTAINS
REAL ddnv(DIMENSION%nv2d,DIMENSION%jspd),dudz(DIMENSION%nv2d,DIMENSION%jspd)
REAL duz(DIMENSION%nv2d,DIMENSION%jspd), udz(DIMENSION%nv2d,DIMENSION%jspd)
REAL uz(DIMENSION%nv2d,DIMENSION%jspd)
! l_J auxiliary potential array
COMPLEX, ALLOCATABLE :: vxy1(:,:,:)
! ..
d2 = SQRT(cell%omtil/cell%area)
IF (jij%l_J) ALLOCATE (vxy1(vacuum%nmzxyd,stars%ng2-1,2))
!---> set up mapping function from 3d-->2d lapws
......
......@@ -9,7 +9,7 @@ MODULE m_od_hsvac
CONTAINS
SUBROUTINE od_hsvac(&
vacuum,stars,DIMENSION, oneD,atoms, jsp,input,vxy,vz,evac,cell,&
bkpt,lapw, MM,vM,m_cyl,n2d_1, n_size,n_rank,sym,noco,jij,nv2,l_real,hamOvlp)
bkpt,lapw, MM,vM,m_cyl,n2d_1, n_size,n_rank,sym,noco,nv2,l_real,hamOvlp)
! subroutine for calculating the hamiltonian and overlap matrices in
! the vacuum in the case of 1-dimensional calculations
......@@ -25,7 +25,6 @@ CONTAINS
TYPE(t_input),INTENT(IN) :: input
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_jij),INTENT(IN) :: jij
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
......@@ -70,14 +69,11 @@ CONTAINS
REAL, ALLOCATABLE :: ddnv(:,:,:),dudz(:,:,:)
REAL, ALLOCATABLE :: duz(:,:,:)
REAL, ALLOCATABLE :: udz(:,:,:),uz(:,:,:)
! l_J auxiliary potential array
COMPLEX, ALLOCATABLE :: vxy1(:,:,:)
! ..
ic = CMPLX(0.,1.)
d2 = SQRT(cell%omtil/cell%area)
IF (jij%l_J) ALLOCATE (vxy1(vacuum%nmzxyd,n2d_1-1,2))
ALLOCATE (&
ai(-vM:vM,DIMENSION%nv2d,DIMENSION%nvd),bi(-vM:vM,DIMENSION%nv2d,DIMENSION%nvd),&
nvp(DIMENSION%nv2d,DIMENSION%jspd),ind(stars%ng2,DIMENSION%nv2d,DIMENSION%jspd),&
......@@ -133,35 +129,13 @@ CONTAINS
!---> load the non-warping part of the potential
READ (25)((vz(imz,ivac,ipot),imz=1,vacuum%nmzd),ipot=1,4)
npot = 3
!---> for J-coeff. we average the up-up and down-down parts
!---> and off-diagonal elements of the potential matrix to zero
IF (jij%l_J) THEN
vz(:,ivac,1) = (vz(:,ivac,1) + vz(:,ivac,2))/2.
vz(:,ivac,2) = vz(:,ivac,1)
vz(:,ivac,3) = 0.0
vz(:,ivac,4) = 0.0
END IF
ENDIF
DO ipot = 1,npot
IF (noco%l_noco) THEN
IF (.NOT.jij%l_J) THEN
READ (25)((vxy(imz,k,ivac), imz=1,vacuum%nmzxy),k=1,n2d_1-1)
END IF
READ (25)((vxy(imz,k,ivac), imz=1,vacuum%nmzxy),k=1,n2d_1-1)
!---> l_J we want to average the diagonal elements of the pot. matrix
IF (jij%l_J .AND. ipot.EQ.1) THEN
READ (25)((vxy(imz,k,ivac), imz=1,vacuum%nmzxy),k=1,n2d_1-1)
READ (25)((vxy1(imz,k,ivac), imz=1,vacuum%nmzxy),k=1,n2d_1-1)
vxy(:,:,ivac) = (vxy(:,:,ivac)+vxy1(:,:,ivac))/2.
END IF
IF (jij%l_J .AND. ipot.EQ.3) THEN
READ (25)((vxy(imz,k,ivac), imz=1,vacuum%nmzxy),k=1,n2d_1-1)
END IF
IF (jij%l_J .AND. ipot.EQ.3) vxy(:,:,ivac)=CMPLX(0.,0.)
ENDIF ! loco
! get the wavefunctions and set up the tuuv, etc matrices
......@@ -354,8 +328,6 @@ CONTAINS
ENDDO !ipot
IF (jij%l_J) DEALLOCATE (vxy1)
DEALLOCATE (ai,bi,nvp,ind,kvac3,map1, tddv,tduv,tudv,tuuv,a,b,bess,dbss,bess1, ddnv,dudz,duz,udz,uz )
RETURN
......
......@@ -11,4 +11,5 @@ eigen_soc/sointg.f90
eigen_soc/sorad.f90
eigen_soc/spnorb.f90
eigen_soc/vso.f90
eigen_soc/ssomat.F90
)
......@@ -13,9 +13,6 @@ MODULE m_spnorb
!*********************************************************************
CONTAINS
SUBROUTINE spnorb(atoms,noco,input,mpi, enpara, vr, usdus, rsoc,l_angles)
USE m_anglso
USE m_sgml
USE m_sorad
USE m_types
IMPLICIT NONE
......@@ -37,12 +34,7 @@ CONTAINS
INTEGER is1,is2,jspin1,jspin2,l,l1,l2,m1,m2,n
LOGICAL, SAVE :: first_k = .TRUE.
! ..
! .. Local Arrays ..
INTEGER ispjsp(2)
! ..
! ..
DATA ispjsp/1,-1/
!Allocate space for SOC matrix elements; set to zero at the same time
ALLOCATE(rsoc%rsopp (atoms%ntype,atoms%lmaxd,2,2));rsoc%rsopp =0.0
ALLOCATE(rsoc%rsoppd (atoms%ntype,atoms%lmaxd,2,2));rsoc%rsoppd=0.0
......@@ -97,10 +89,34 @@ CONTAINS
9000 FORMAT (5x,' p ',5x,' d ', 5x, ' f ')
!
IF (.NOT.l_angles) RETURN
!Calculate angular matrix elements if requested
IF (l_angles) &
CALL spnorb_angles(atoms,mpi,noco%theta,noco%phi,rsoc%soangl)
END SUBROUTINE spnorb
IF ((ABS(noco%theta).LT.0.00001).AND.(ABS(noco%phi).LT.0.00001)) THEN
SUBROUTINE spnorb_angles(atoms,mpi,theta,phi,soangl)
USE m_anglso
USE m_sgml
USE m_sorad
USE m_types
IMPLICIT NONE
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_mpi),INTENT(IN) :: mpi
REAL,INTENT(IN) :: theta,phi
COMPLEX,INTENT(INOUT) :: soangl(:,-atoms%lmaxd:,:,:,-atoms%lmaxd:,:)
! ..
! ..
! .. Local Scalars ..
INTEGER is1,is2,jspin1,jspin2,l,l1,l2,m1,m2,n
! ..
! .. Local Arrays ..
INTEGER ispjsp(2)
! ..
! ..
DATA ispjsp/1,-1/
IF ((ABS(theta).LT.0.00001).AND.(ABS(phi).LT.0.00001)) THEN
!
! TEST for real function sgml(l1,m1,is1,l2,m2,is2)
!
......@@ -112,7 +128,7 @@ CONTAINS
is2=ispjsp(jspin2)
DO m1 = -l1,l1,1
DO m2 = -l2,l2,1
rsoc%soangl(l1,m1,jspin1,l2,m2,jspin2) =&
soangl(l1,m1,jspin1,l2,m2,jspin2) =&
CMPLX(sgml(l1,m1,is1,l2,m2,is2),0.0)
ENDDO
ENDDO
......@@ -134,8 +150,8 @@ CONTAINS
!
DO m1 = -l1,l1,1
DO m2 = -l2,l2,1
rsoc%soangl(l1,m1,jspin1,l2,m2,jspin2) =&
anglso(noco%theta,noco%phi,l1,m1,is1,l2,m2,is2)
soangl(l1,m1,jspin1,l2,m2,jspin2) =&
anglso(theta,phi,l1,m1,is1,l2,m2,is2)
ENDDO
ENDDO
!
......@@ -152,7 +168,7 @@ CONTAINS
DO jspin2 = 1,2
WRITE (6,FMT=*) 'd-states:is1=',jspin1,',is2=',jspin2
WRITE (6,FMT='(7x,7i8)') (m1,m1=-3,3,1)
WRITE (6,FMT=8003) (m2, (rsoc%soangl(3,m1,jspin1,3,m2,jspin2),&
WRITE (6,FMT=8003) (m2, (soangl(3,m1,jspin1,3,m2,jspin2),&
m1=-3,3,1),m2=-3,3,1)
ENDDO
ENDDO
......@@ -160,5 +176,5 @@ CONTAINS
8002 FORMAT (' so - angular matrix elements')
8003 FORMAT (i8,14f8.4)
END SUBROUTINE spnorb
END SUBROUTINE spnorb_angles
END MODULE m_spnorb
MODULE m_ssomat
USE m_judft
IMPLICIT NONE
CONTAINS
SUBROUTINE ssomat(seigvso,theta,phi,eig_id,DIMENSION,atoms,kpts,sym,&
cell,noco, input,mpi, oneD,enpara,v,results )
USE m_eig66_io
USE m_spnorb
USE m_abcof
USE m_types_setup
USE m_types_mpi
USE m_types_enpara
USE m_types_potden
USE m_types_misc
USE m_types_kpts
USE m_types_tlmplm
USE m_types_usdus
USE m_types_lapw
USE m_constants
IMPLICIT NONE
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_dimension),INTENT(IN) :: dimension
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_enpara),INTENT(IN) :: enpara
TYPE(t_potden),INTENT(IN) :: v
TYPE(t_results),INTENT(IN) :: results
INTEGER,INTENT(IN) :: eig_id
REAL,INTENT(in) :: theta(:),phi(:) ! more than a single angle at once...
REAL,INTENT(OUT) :: seigvso(:)
! ..
! .. Locals ..
#ifdef CPP_MPI
INTEGER:: ierr
include 'mpif.h'
#endif
INTEGER :: neigf=1 !not full-matrix
INTEGER :: ilo,js,jsloc,nk,n,l ,lm,band,nr,ne,nat,m
INTEGER :: na
REAL :: r1,r2
COMPLEX :: c1,c2
COMPLEX, ALLOCATABLE :: matel(:,:,:)
REAL, ALLOCATABLE :: eig_shift(:,:,:)
COMPLEX, ALLOCATABLE :: acof(:,:,:,:,:), bcof(:,:,:,:,:)
COMPLEX, ALLOCATABLE :: ccof(:,:,:,:,:,:)
COMPLEX,ALLOCATABLE :: soangl(:,:,:,:,:,:,:)
TYPE(t_rsoc):: rsoc
TYPE(t_zmat):: zmat
TYPE(t_noco):: noco_tmp
TYPE(t_usdus):: usdus
TYPE(t_lapw) :: lapw
IF (ANY(atoms%neq/=1)) CALL judft_error('(spin spiral + soc) does not work'//&
' properly for more than one atom per type!',calledby="ssomat")
! needed directly for calculating matrix elements
seigvso=0.0
ALLOCATE(eig_shift(DIMENSION%neigd,kpts%nkpt,SIZE(theta)));eig_shift=0.0
ALLOCATE( acof(dimension%neigd,0:dimension%lmd,atoms%nat,2,2),&
bcof(dimension%neigd,0:dimension%lmd,atoms%nat,2,2) )
ALLOCATE( ccof(-atoms%llod:atoms%llod,dimension%neigd,atoms%nlod,atoms%nat,2,2) )
ALLOCATE( matel(neigf,DIMENSION%neigd,0:atoms%ntype) )
zMat%nbasfcn=lapw%nv(1)+atoms%nlotot
zmat%nbands=DIMENSION%neigd
zmat%l_real=.FALSE.
ALLOCATE(zmat%z_c(zMat%nbasfcn,zmat%nbands))
CALL usdus%init(atoms,2)
!Calculate radial and angular matrix elements of SOC
!many directions of SOC at once...
CALL spnorb(atoms,noco_tmp,input,mpi, enpara, v%mt, usdus, rsoc,.FALSE.)
ALLOCATE(soangl(atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,2,&
atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,2,SIZE(theta)))
DO nr=1,SIZE(theta)
CALL spnorb_angles(atoms,mpi,theta(nr),phi(nr),soangl(:,:,:,:,:,:,nr))
ENDDO
DO nk=mpi%irank+1,kpts%nkpt,mpi%isize
CALL lapw%init(input,noco, kpts,atoms,sym,nk,cell,.false., mpi)
DO jsloc= 1,2
CALL read_eig(eig_id,nk,jsloc,neig=ne,zmat=zmat)
CALL abcof(input,atoms,sym, cell,lapw,ne,usdus,noco,jsloc,oneD, &
acof(:,:,:,jsloc,1),bcof(:,:,:,jsloc,1),ccof(:,:,:,:,jsloc,1),zMat)
ENDDO
! rotate abcof into global spin coordinate frame
nat= 0
DO n= 1,atoms%ntype
DO na= 1,atoms%neq(n)
nat= nat+1
r1= noco%alph(n)
r2= noco%beta(n)
DO lm= 0,DIMENSION%lmd
DO band= 1,DIMENSION%neigd
c1= acof(band,lm,nat,1,1)
c2= acof(band,lm,nat,2,1)
acof(band,lm,nat,1,1)= CMPLX(COS(r1/2.),-SIN(r1/2.)) *CMPLX( COS(r2/2.),0.) *c1
acof(band,lm,nat,2,1)= CMPLX(COS(r1/2.),-SIN(r1/2.)) *CMPLX(-SIN(r2/2.),0.) *c2
acof(band,lm,nat,1,2)= CMPLX(COS(r1/2.),+SIN(r1/2.)) *CMPLX(+SIN(r2/2.),0.) *c1
acof(band,lm,nat,2,2)= CMPLX(COS(r1/2.),+SIN(r1/2.)) *CMPLX( COS(r2/2.),0.) *c2
c1= bcof(band,lm,nat,1,1)
c2= bcof(band,lm,nat,2,1)
bcof(band,lm,nat,1,1)= CMPLX(COS(r1/2.),-SIN(r1/2.)) *CMPLX( COS(r2/2.),0.) *c1
bcof(band,lm,nat,2,1)= CMPLX(COS(r1/2.),-SIN(r1/2.)) *CMPLX(-SIN(r2/2.),0.) *c2
bcof(band,lm,nat,1,2)= CMPLX(COS(r1/2.),+SIN(r1/2.)) *CMPLX(+SIN(r2/2.),0.) *c1
bcof(band,lm,nat,2,2)= CMPLX(COS(r1/2.),+SIN(r1/2.)) *CMPLX( COS(r2/2.),0.) *c2
ENDDO ! band
ENDDO ! lm
DO ilo = 1,atoms%nlo(n)
l = atoms%llo(ilo,n)
DO band= 1,DIMENSION%neigd
DO m = -l, l
c1= ccof(m,band,ilo,nat,1,1)
c2= ccof(m,band,ilo,nat,2,1)
ccof(m,band,ilo,nat,1,1)= CMPLX(COS(r1/2.),-SIN(r1/2.))*CMPLX( COS(r2/2.),0.)*c1
ccof(m,band,ilo,nat,2,1)= CMPLX(COS(r1/2.),-SIN(r1/2.))*CMPLX(-SIN(r2/2.),0.)*c2
ccof(m,band,ilo,nat,1,2)= CMPLX(COS(r1/2.),+SIN(r1/2.))*CMPLX(+SIN(r2/2.),0.)*c1
ccof(m,band,ilo,nat,2,2)= CMPLX(COS(r1/2.),+SIN(r1/2.))*CMPLX( COS(r2/2.),0.)*c2
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
DO nr=1,size(theta) !loop over angles
! matrix elements within k
CALL ssomatel(neigf,dimension,atoms, noco, &
soangl(:,:,:,:,:,:,nr),rsoc%rsopp(:,:,:,:),rsoc%rsoppd(:,:,:,:),&
rsoc%rsopdp(:,:,:,:),rsoc%rsopdpd(:,:,:,:),rsoc%rsoplop(:,:,:,:), &
rsoc%rsoplopd(:,:,:,:),rsoc%rsopdplo(:,:,:,:),rsoc%rsopplo(:,:,:,:),&
rsoc%rsoploplop(:,:,:,:,:),&
.TRUE.,&
acof,bcof, ccof,&
acof,bcof, ccof,&
matel )
eig_shift(:,nk,nr)=matel(1,:,0)
ENDDO
ENDDO
!Collect data from distributed k-loop
#ifdef CPP_MPI
IF (mpi%irank==0) THEN
CALL MPI_REDUCE(MPI_IN_PLACE,eig_shift,SIZE(eig_shift),MPI_DOUBLE_PRECISION,MPI_SUM,0,mpi%mpi_comm,ierr)
ELSE
CALL MPI_REDUCE(eig_shift,eig_shift,SIZE(eig_shift),MPI_DOUBLE_PRECISION,MPI_SUM,0,mpi%mpi_comm,ierr)
ENDIF
#endif
IF (mpi%irank==0) THEN
!Sum all shift using weights
DO nr=1,SIZE(theta)
DO nk=1,kpts%nkpt
seigvso(nr)=seigvso(nr)+dot_PRODUCT(results%w_iks(:,nk,1),eig_shift(:,nk,nr))
ENDDO
ENDDO
seigvso= results%seigv+seigvso
ENDIF
END SUBROUTINE ssomat
! ==================================================================== !
SUBROUTINE ssomatel(neigf,dimension,atoms, noco,&
soangl,rsopp,rsoppd,rsopdp,rsopdpd,rsoplop,&
rsoplopd,rsopdplo,rsopplo,rsoploplop,&
diag, &
acof1,bcof1,ccof1,acof2,bcof2,ccof2,&
matel )
USE m_types
IMPLICIT NONE
TYPE(t_dimension),INTENT(IN) :: dimension
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_atoms),INTENT(IN) :: atoms
LOGICAL, INTENT(IN) :: diag
INTEGER, INTENT(IN) :: neigf
REAL, INTENT(IN) :: &
rsopp(:,:,:,:), rsoppd(:,:,:,:),&
rsopdp(:,:,:,:), rsopdpd(:,:,:,:), &
rsoplop(:,:,:,:),rsoplopd(:,:,:,:),&
rsopdplo(:,:,:,:),rsopplo(:,:,:,:),&
rsoploplop(:,:,:,:,:)
COMPLEX, INTENT(IN) :: &
soangl(:,-atoms%lmaxd:,:,:,-atoms%lmaxd:,:), &
acof1(:,0:,:,:,:), &
bcof1(:,0:,:,:,:),&
ccof1(-atoms%llod:,:,:,:,:,:),&
acof2(:,0:,:,:,:), &
bcof2(:,0:,:,:,:),&
ccof2(-atoms%llod:,:,:,:,:,:)
Complex, INTENT(OUT) :: matel(neigf,dimension%neigd,0:atoms%ntype)
INTEGER :: band1,band2,bandf, n ,na, l,m1,m2,lm1,lm2,&
jsloc1,jsloc2, js1,js2,jsnumber,ilo,ilop,nat
COMPLEX, ALLOCATABLE :: sa(:,:),sb(:,:),sc(:,:,:),ral(:,:,:)
COMPLEX, ALLOCATABLE :: ra(:,:),rb(:,:),rc(:,:,:),rbl(:,:,:)
! with the following nesting of loops the calculation of the
! matrix-elements is of order
! natall*lmd*neigd*(lmd+neigd) ; note that lmd+neigd << lmd*neigd
matel(:,:,:)= CMPLX(0.,0.)
ALLOCATE ( sa(2,0:dimension%lmd),sb(2,0:dimension%lmd),ra(2,0:dimension%lmd),rb(2,0:dimension%lmd) )
ALLOCATE ( sc(2,-atoms%llod:atoms%llod,atoms%nlod),rc(2,-atoms%llod:atoms%llod,atoms%nlod) )
ALLOCATE ( ral(2,-atoms%llod:atoms%llod,atoms%nlod),rbl(2,-atoms%llod:atoms%llod,atoms%nlod) )
! within one k-point loop over global spin
IF (diag) THEN
jsnumber= 2
ELSE
jsnumber= 1
ENDIF
DO js2= 1,jsnumber
IF (diag) THEN
js1= js2
ELSE
js1= 2
ENDIF
! loop over MT
na= 0
DO n= 1,atoms%ntype
DO nat= 1,atoms%neq(n)
na= na+1
DO band2= 1,dimension%neigd ! loop over eigenstates 2
DO l= 1,atoms%lmax(n) ! loop over l
DO m1= -l,l ! loop over m1
lm1= l*(l+1) + m1
DO jsloc2= 1,2
sa(jsloc2,lm1) = CMPLX(0.,0.)
sb(jsloc2,lm1) = CMPLX(0.,0.)
DO m2= -l,l
lm2= l*(l+1) + m2
sa(jsloc2,lm1)= sa(jsloc2,lm1) + &
CONJG(acof2(band2,lm2,na,jsloc2,js2))&
* soangl(l,m2,js2,l,m1,js1)
sb(jsloc2,lm1)= sb(jsloc2,lm1) + &
CONJG(bcof2(band2,lm2,na,jsloc2,js2))&
* soangl(l,m2,js2,l,m1,js1)
ENDDO ! m2
ENDDO ! jsloc2
ENDDO ! m1
ENDDO ! l
DO ilo = 1, atoms%nlo(n) ! LO-part
l = atoms%llo(ilo,n)
DO m1 = -l, l
DO jsloc2= 1,2
sc(jsloc2,m1,ilo) = CMPLX(0.,0.)
DO m2= -l, l
sc(jsloc2,m1,ilo) = sc(jsloc2,m1,ilo) +&
CONJG(ccof2(m2,band2,ilo,na,jsloc2,js2))&
* soangl(l,m2,js2,l,m1,js1)
ENDDO
ENDDO
ENDDO
ENDDO ! ilo
DO l= 1,atoms%lmax(n) ! loop over l
DO m1= -l,l ! loop over m1
lm1= l*(l+1) + m1
DO jsloc1= 1,2
ra(jsloc1,lm1)= CMPLX(0.,0.)
rb(jsloc1,lm1)= CMPLX(0.,0.)
DO jsloc2= 1,2