Commit 59e03acc authored by Daniel Wortmann's avatar Daniel Wortmann

Further refactoring of Wannier code. Should compile and run simple example~

parent 13992254
......@@ -4,9 +4,7 @@ wannier/eulerrot.f
#wannier/w90kpunktgen.f
wannier/wann_1dvacabcof.F
wannier/wann_2dvacabcof.F
wannier/wann_abinv.f
wannier/wann_amn.f
wannier/wann_anglmom.f
wannier/wann_dipole2.f
wannier/wann_dipole3.f
wannier/wann_dipole_electronic.f
......@@ -27,7 +25,6 @@ wannier/wann_gwf_commat.f
wannier/wann_gwf_tools.f
wannier/wann_gwf_write_mmnk.F
wannier/wann_hopping.F
wannier/wannier.F
wannier/wannier_to_lapw.F
wannier/wann_ioncharge_gen.f
wannier/wann_kpointgen.f
......@@ -52,7 +49,6 @@ wannier/wann_mmnk_symm.f
wannier/wann_nabla_pauli_rs.f
wannier/wann_nabla_rs.f
wannier/wann_nocoplot.F
wannier/wann_orbcomp.f
#wannier/wann_orbmag.F
wannier/wann_pauli_rs.F
wannier/wann_perpmag_rs.f
......@@ -62,7 +58,6 @@ wannier/wann_plot_symm.f
wannier/wann_plot_um_dat.F
wannier/wann_plot_vac.F
#wannier/wann_plotw90.F
wannier/wann_postproc.F
wannier/wann_postproc_setup4.F
wannier/wann_postproc_setup5.F
wannier/wann_postproc_setup.F
......@@ -99,6 +94,11 @@ set(fleur_F90 ${fleur_F90}
wannier/init_wannier_defaults.f90
wannier/wann_read_inp.f90
wannier/wann_optional.f90
wannier/wann_abinv.f90
wannier/wann_anglmom.f90
wannier/wannier.F90
wannier/wann_orbcomp.f90
wannier/wann_postproc.F90
)
if(FLEUR_USE_WANN)
set(fleur_F90 ${fleur_F90}
......
......@@ -1100,9 +1100,7 @@ c***********************************************************
DEALLOCATE(lapw_b%k1,lapw_b%k2,lapw_b%k3)
call wann_abinv(
> ntypd,natd,noccbd_b,lmaxd,lmd,llod,nlod,ntype,neq,
> noccbd_b,lmax,nlo,llo,invsat,invsatnr,bkpt_b,taual,
call wann_abinv(atoms,
X acof_b,bcof_b,ccof_b)
call cpu_time(t1)
t_abcof = t_abcof + t1 - t0
......@@ -1222,9 +1220,7 @@ c***********************************************************
DEALLOCATE(lapw_b2%k1,lapw_b2%k2,lapw_b2%k3)
call wann_abinv(
> ntypd,natd,noccbd_b2,lmaxd,lmd,llod,nlod,ntype,neq,
> noccbd_b2,lmax,nlo,llo,invsat,invsatnr,bkpt_b2,taual,
call wann_abinv(atoms,
X acof_b2,bcof_b2,ccof_b2)
call cpu_time(t1)
t_abcof = t_abcof + t1 - t0
......
......@@ -1073,9 +1073,7 @@ c***********************************************************
DEALLOCATE(lapw_b%k1,lapw_b%k2,lapw_b%k3)
call wann_abinv(
> ntypd,natd,noccbd_b,lmaxd,lmd,llod,nlod,ntype,neq,
> noccbd_b,lmax,nlo,llo,invsat,invsatnr,bkpt_b,taual,
call wann_abinv(atoms,
X acof_b,bcof_b,ccof_b)
call cpu_time(t1)
t_abcof = t_abcof + t1 - t0
......@@ -1205,9 +1203,7 @@ c***********************************************************
DEALLOCATE(lapw_b2%k1,lapw_b2%k2,lapw_b2%k3)
call wann_abinv(
> ntypd,natd,noccbd_b2,lmaxd,lmd,llod,nlod,ntype,neq,
> noccbd_b2,lmax,nlo,llo,invsat,invsatnr,bkpt_b2,taual,
call wann_abinv(atoms,
X acof_b2,bcof_b2,ccof_b2)
call cpu_time(t1)
t_abcof = t_abcof + t1 - t0
......
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
module m_wann_abinv
contains
SUBROUTINE wann_abinv(
> ntypd,natd,neigd,lmaxd,lmd,llod,nlod,ntype,neq,
> neig,lmax,nlo,llo,invsat,invsatnr,bkpt,taual,
X acof,bcof,ccof)
C ***************************************************************
C Transform acof,bcof,ccof in case of atoms related by inversion
c symmetry to obtain the coefficients in the global frame.
c Based on abcrot.
c Frank Freimuth
C ***************************************************************
use m_constants, only:pimach
IMPLICIT NONE
C ..
C .. Scalar Arguments ..
INTEGER, INTENT (IN) :: ntypd,natd,neigd,lmd,llod,nlod,ntype
INTEGER, INTENT (IN) :: lmaxd,neig
C ..
C .. Array Arguments ..
INTEGER, INTENT (IN) :: neq(ntypd),lmax(ntypd),nlo(ntypd)
INTEGER, INTENT (IN) :: llo(nlod,ntypd)
INTEGER, INTENT (IN) :: invsat(natd),invsatnr(natd)
real,intent(in) :: bkpt(3)
REAL, INTENT (IN) :: taual(3,natd)
COMPLEX, INTENT (INOUT) :: acof(neigd,0:lmd,natd)
COMPLEX, INTENT (INOUT) :: bcof(neigd,0:lmd,natd)
COMPLEX, INTENT (INOUT) :: ccof(-llod:llod,neigd,nlod,natd)
C .. Local Scalars ..
INTEGER :: itype,ineq,iatom,iop,ilo,i,l,m,lm,lmp,ifac
integer :: n,nn,jatom,ie,ll1
real :: tpi,arg
complex :: fac
C ..
tpi=2.0*pimach()
iatom=0
DO itype=1,ntype
DO ineq=1,neq(itype)
iatom=iatom+1
IF(invsat(iatom).ne.2) cycle
DO l=1,lmax(itype),2
DO i=1,neig
acof(i,l**2:l*(l+2),iatom) = (-1)**l *
& acof(i,l**2:l*(l+2),iatom)
bcof(i,l**2:l*(l+2),iatom) = (-1)**l *
& bcof(i,l**2:l*(l+2),iatom)
ENDDO
ENDDO
DO ilo=1,nlo(itype)
l=llo(ilo,itype)
IF(l.gt.0) THEN
if(mod(l,2).eq.0)cycle
DO i=1,neig
ccof(-l:l,i,ilo,iatom) = (-1)**l *
& ccof(-l:l,i,ilo,iatom)
ENDDO
ENDIF
ENDDO
ENDDO
ENDDO
c$$$ iatom = 0
c$$$ DO n = 1,ntype
c$$$ DO nn = 1,neq(n)
c$$$ iatom = iatom + 1
c$$$ IF (invsat(iatom).EQ.1) THEN
c$$$ jatom = invsatnr(iatom)
c$$$ arg= (taual(1,jatom)+taual(1,iatom))*bkpt(1)
c$$$ arg=arg+(taual(2,jatom)+taual(2,iatom))*bkpt(2)
c$$$ arg=arg+(taual(3,jatom)+taual(3,iatom))*bkpt(3)
c$$$ arg=arg*tpi
c$$$ fac=cmplx(cos(arg),sin(arg))
c$$$ DO ilo = 1,nlo(n)
c$$$ l = llo(ilo,n)
c$$$ DO m = -l,l
c$$$ DO ie = 1,neig
c$$$ ccof(m,ie,ilo,jatom) = fac *
c$$$ + ccof(m,ie,ilo,jatom)
c$$$ ENDDO
c$$$ ENDDO
c$$$ ENDDO
c$$$ DO l = 0,lmax(n)
c$$$ ll1 = l* (l+1)
c$$$ DO m =-l,l
c$$$ lm = ll1 + m
c$$$ DO ie = 1,neig
c$$$ acof(ie,lm,jatom) = fac *
c$$$ * acof(ie,lm,jatom)
c$$$ ENDDO
c$$$ DO ie = 1,neig
c$$$ bcof(ie,lm,jatom) = fac *
c$$$ * bcof(ie,lm,jatom)
c$$$ ENDDO
c$$$ ENDDO
c$$$ ENDDO
c$$$ ENDIF
c$$$ ENDDO
c$$$ ENDDO
END subroutine
end module
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_wann_abinv
CONTAINS
SUBROUTINE wann_abinv(atoms,acof,bcof,ccof)
! ***************************************************************
! Transform acof,bcof,ccof in case of atoms related by inversion
! symmetry to obtain the coefficients in the global frame.
! Based on abcrot.
! Frank Freimuth
! ***************************************************************
USE m_types
IMPLICIT NONE
! ..
! .. Scalar Arguments ..
TYPE(t_atoms),INTENT(IN) :: atoms
COMPLEX, INTENT (INOUT) :: acof(:,0:,:)
COMPLEX, INTENT (INOUT) :: bcof(:,0:,:)
COMPLEX, INTENT (INOUT) :: ccof(-atoms%llod:,:,:,:)!(-llod:llod,neigd,nlod,natd)
! .. Local Scalars ..
INTEGER :: itype,ineq,iatom,ilo,l
iatom=0
DO itype=1,atoms%ntype
DO ineq=1,atoms%neq(itype)
iatom=iatom+1
IF(atoms%invsat(iatom).NE.2) CYCLE
DO l=1,atoms%lmax(itype),2
acof(:,l**2:l*(l+2),iatom) = (-1)**l *&
acof(:,l**2:l*(l+2),iatom)
bcof(:,l**2:l*(l+2),iatom) = (-1)**l * &
bcof(:,l**2:l*(l+2),iatom)
ENDDO
DO ilo=1,atoms%nlo(itype)
l=atoms%llo(ilo,itype)
IF(l.GT.0) THEN
IF(MOD(l,2).EQ.0)CYCLE
ccof(-l:l,:,ilo,iatom) = (-1)**l * &
ccof(-l:l,:,ilo,iatom)
ENDIF
ENDDO
ENDDO
ENDDO
END SUBROUTINE wann_abinv
END MODULE m_wann_abinv
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_wann_anglmom
c***********************************************************************
c Compute matrix elements of angular momentum operator
c in the muffin-tin spheres.
c
c Frank Freimuth
c***********************************************************************
CONTAINS
SUBROUTINE wann_anglmom(
> llod,noccbd,nlod,natd,ntypd,lmax,lmd,
> ntype,neq,nlo,llo,acof,bcof,ccof,
> ddn,uulon,dulon,uloulopn,
= mmn)
implicit none
c .. scalar arguments ..
integer, intent (in) :: llod,nlod,natd,ntypd,lmd
integer, intent (in) :: ntype,noccbd
c .. array arguments ..
integer, intent (in) :: neq(:)!neq(ntypd)
integer, intent (in) :: nlo(:)!nlo(ntypd)
integer, intent (in) :: lmax(:)!lmax(ntypd)
integer, intent (in) :: llo(:,:)!llo(nlod,ntypd)
real, intent (in) :: ddn(0:,:)!ddn(0:lmaxd,ntypd)
real, intent (in) :: uloulopn(:,:,:)!uloulopn(nlod,nlod,ntypd)
real, intent (in) :: uulon(:,:)!uulon(nlod,ntypd)
real, intent (in) :: dulon(:,:)!dulon(nlod,ntypd)
complex, intent (in) :: ccof(-llod:,:,:,:) !ccof(-llod:llod,noccbd,nlod,natd)
complex, intent (in) :: acof(:,0:,:)!acof(noccbd,0:lmd,natd)
complex, intent (in) :: bcof(:,0:,:)!bcof(noccbd,0:lmd,natd)
complex, intent (inout) :: mmn(:,:,:)!mmn(3,noccbd,noccbd)
c .. local scalars ..
logical :: l_select
integer :: i,j,l,lo,lop,m,natom,nn,ntyp
integer :: nt1,nt2,lm,n,ll1,indat
complex :: suma_z,sumb_z
complex :: suma_p,sumb_p
complex :: suma_m,sumb_m
complex :: suma_x,sumb_x
complex :: suma_y,sumb_y
real :: lplus,lminus
C ..
C .. local arrays ..
complex, allocatable :: qlo_z(:,:,:,:,:)
complex, allocatable :: qlo_p(:,:,:,:,:)
complex, allocatable :: qlo_m(:,:,:,:,:)
complex, allocatable :: qaclo_z(:,:,:,:),qbclo_z(:,:,:,:)
complex, allocatable :: qaclo_p(:,:,:,:),qbclo_p(:,:,:,:)
complex, allocatable :: qaclo_m(:,:,:,:),qbclo_m(:,:,:,:)
C ..
C .. intrinsic functions ..
intrinsic conjg
allocate (qlo_z(noccbd,noccbd,nlod,nlod,ntypd),
+ qaclo_z(noccbd,noccbd,nlod,ntypd),
+ qbclo_z(noccbd,noccbd,nlod,ntypd) )
allocate (qlo_p(noccbd,noccbd,nlod,nlod,ntypd),
+ qaclo_p(noccbd,noccbd,nlod,ntypd),
+ qbclo_p(noccbd,noccbd,nlod,ntypd) )
allocate (qlo_m(noccbd,noccbd,nlod,nlod,ntypd),
+ qaclo_m(noccbd,noccbd,nlod,ntypd),
+ qbclo_m(noccbd,noccbd,nlod,ntypd) )
inquire(file='select_anglmom',exist=l_select)
write(*,*)'select_anglmom: ',l_select
if(l_select) then
open(866,file='select_anglmom')
read(866,*)indat
close(866)
write(*,*)'anglmom for atom=',indat
write(*,*)ntype
write(*,*)neq(indat)
endif
c-----> lapw-lapw-Terms
do i = 1,noccbd
do j = 1,noccbd
nt1 = 1
do n = 1,ntype
nt2 = nt1 + neq(n) - 1
do l = 0,lmax(n)
suma_z = cmplx(0.,0.); sumb_z = cmplx(0.,0.)
suma_m = cmplx(0.,0.); sumb_m = cmplx(0.,0.)
suma_p = cmplx(0.,0.); sumb_p = cmplx(0.,0.)
if(l_select .and. (n.ne.indat)) cycle
ll1 = l* (l+1)
do m = -l,l
lm = ll1 + m
lplus=sqrt(real( (l-m)*(l+m+1) ) )
lminus=sqrt(real( (l+m)*(l-m+1) ) )
do natom = nt1,nt2
suma_z = suma_z + acof(i,lm,natom)*
+ conjg(acof(j,lm,natom))*real(m)
sumb_z = sumb_z + bcof(i,lm,natom)*
+ conjg(bcof(j,lm,natom))*real(m)
if(m+1.le.l)then
suma_p = suma_p + acof(i,lm,natom)*
+ conjg(acof(j,lm+1,natom))*lplus
sumb_p = sumb_p + bcof(i,lm,natom)*
+ conjg(bcof(j,lm+1,natom))*lplus
endif
if(m-1.ge.-l)then
suma_m = suma_m + acof(i,lm,natom)*
+ conjg(acof(j,lm-1,natom))*lminus
sumb_m = sumb_m + bcof(i,lm,natom)*
+ conjg(bcof(j,lm-1,natom))*lminus
endif
enddo
enddo
mmn(3,j,i) = mmn(3,j,i) + (suma_z+sumb_z*ddn(l,n))
suma_x=0.5*(suma_p+suma_m)
sumb_x=0.5*(sumb_p+sumb_m)
mmn(1,j,i) = mmn(1,j,i) + (suma_x+sumb_x*ddn(l,n))
suma_y=cmplx(0.0,-0.5)*(suma_p-suma_m)
sumb_y=cmplx(0.0,-0.5)*(sumb_p-sumb_m)
mmn(2,j,i) = mmn(2,j,i) + (suma_y+sumb_y*ddn(l,n))
enddo ! l
nt1 = nt1 + neq(n)
enddo ! n
enddo ! j
enddo ! i
c---> Terms involving local orbitals.
qlo_z = 0.0; qlo_p = 0.0; qlo_m = 0.0
qaclo_z = 0.0; qaclo_p = 0.0; qaclo_m = 0.0
qbclo_z = 0.0; qbclo_p = 0.0; qbclo_m = 0.0
natom = 0
do ntyp = 1,ntype
do nn = 1,neq(ntyp)
natom = natom + 1
if(l_select .and. (ntyp.ne.indat)) cycle
do lo = 1,nlo(ntyp)
l = llo(lo,ntyp)
ll1 = l* (l+1)
do m = -l,l
lm = ll1 + m
lplus=sqrt(real( (l-m)*(l+m+1) ) )
lminus=sqrt(real( (l+m)*(l-m+1) ) )
do i = 1,noccbd
do j = 1,noccbd
qbclo_z(j,i,lo,ntyp) = qbclo_z(j,i,lo,ntyp) + (
+ bcof(i,lm,natom) * conjg(ccof(m,j,lo,natom)) +
+ ccof(m,i,lo,natom)*conjg(bcof(j,lm,natom)) )*real(m)
qaclo_z(j,i,lo,ntyp) = qaclo_z(j,i,lo,ntyp) + (
+ acof(i,lm,natom) * conjg(ccof(m,j,lo,natom)) +
+ ccof(m,i,lo,natom)*conjg(acof(j,lm,natom)) )*real(m)
if(m+1.le.l)then
qbclo_p(j,i,lo,ntyp) = qbclo_p(j,i,lo,ntyp) + (
+ bcof(i,lm,natom) * conjg(ccof(m+1,j,lo,natom)) +
+ ccof(m,i,lo,natom)*conjg(bcof(j,lm+1,natom)) )*lplus
qaclo_p(j,i,lo,ntyp) = qaclo_p(j,i,lo,ntyp) + (
+ acof(i,lm,natom) * conjg(ccof(m+1,j,lo,natom)) +
+ ccof(m,i,lo,natom)*conjg(acof(j,lm+1,natom)) )*lplus
endif
if(m-1.ge.-l)then
qbclo_m(j,i,lo,ntyp) = qbclo_m(j,i,lo,ntyp) + (
+ bcof(i,lm,natom) * conjg(ccof(m-1,j,lo,natom)) +
+ ccof(m,i,lo,natom)*conjg(bcof(j,lm-1,natom)) )*lminus
qaclo_m(j,i,lo,ntyp) = qaclo_m(j,i,lo,ntyp) + (
+ acof(i,lm,natom) * conjg(ccof(m-1,j,lo,natom)) +
+ ccof(m,i,lo,natom)*conjg(acof(j,lm-1,natom)) )*lminus
endif
enddo !j
enddo !i
enddo !m
do lop = 1,nlo(ntyp)
if (llo(lop,ntyp).eq.l) then
do m = -l,l
lplus=sqrt(real( (l-m)*(l+m+1) ) )
lminus=sqrt(real( (l+m)*(l-m+1) ) )
do i = 1,noccbd
do j = 1,noccbd
qlo_z(j,i,lop,lo,ntyp) = qlo_z(j,i,lop,lo,ntyp) +
+ conjg(ccof(m,j,lop,natom))
* *ccof(m,i,lo,natom)*real(m)
if(m+1.le.l)then
qlo_p(j,i,lop,lo,ntyp) =
+ qlo_p(j,i,lop,lo,ntyp) +
+ conjg(ccof(m+1,j,lop,natom))
* *ccof(m,i,lo,natom)*lplus
endif
if(m-1.ge.-l)then
qlo_m(j,i,lop,lo,ntyp) =
+ qlo_m(j,i,lop,lo,ntyp) +
+ conjg(ccof(m-1,j,lop,natom))
* *ccof(m,i,lo,natom)*lminus
endif
enddo ! j
enddo ! i
enddo ! m
endif
enddo ! lop
enddo ! lo
enddo ! nn
enddo ! ntyp
c---> perform summation of the coefficients with the integrals
c---> of the radial basis functions
do ntyp = 1,ntype
if(l_select .and. (ntyp.ne.indat) ) cycle
do lo = 1,nlo(ntyp)
l = llo(lo,ntyp)
do j = 1,noccbd
do i = 1,noccbd
mmn(3,i,j)= mmn(3,i,j) +
+ qaclo_z(i,j,lo,ntyp)*uulon(lo,ntyp) +
+ qbclo_z(i,j,lo,ntyp)*dulon(lo,ntyp)
suma_p=qaclo_p(i,j,lo,ntyp)*uulon(lo,ntyp) +
+ qbclo_p(i,j,lo,ntyp)*dulon(lo,ntyp)
suma_m=qaclo_m(i,j,lo,ntyp)*uulon(lo,ntyp) +
+ qbclo_m(i,j,lo,ntyp)*dulon(lo,ntyp)
suma_x= 0.5*(suma_p+suma_m)
suma_y=cmplx(0.0,-0.5)*(suma_p-suma_m)
mmn(1,i,j)= mmn(1,i,j) + suma_x
mmn(2,i,j)= mmn(2,i,j) + suma_y
enddo !i
enddo !j
do lop = 1,nlo(ntyp)
if (llo(lop,ntyp).eq.l) then
do j = 1,noccbd
do i = 1,noccbd
mmn(3,i,j) = mmn(3,i,j) +
+ qlo_z(i,j,lop,lo,ntyp)*uloulopn(lop,lo,ntyp)
suma_p=qlo_p(i,j,lop,lo,ntyp)*uloulopn(lop,lo,ntyp)
suma_m=qlo_m(i,j,lop,lo,ntyp)*uloulopn(lop,lo,ntyp)
mmn(1,i,j) = mmn(1,i,j) + 0.5*(suma_p+suma_m)
mmn(2,i,j) = mmn(2,i,j) +
+ cmplx(0.0,-0.5)*(suma_p-suma_m)
enddo ! i
enddo ! j
endif
enddo !lop
enddo !lo
enddo !ntyp
deallocate ( qlo_z,qaclo_z,qbclo_z )
deallocate ( qlo_m,qaclo_m,qbclo_m )
deallocate ( qlo_p,qaclo_p,qbclo_p )
END SUBROUTINE wann_anglmom
END MODULE m_wann_anglmom
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_wann_anglmom
!***********************************************************************
! Compute matrix elements of angular momentum operator
! in the muffin-tin spheres.
!
! Frank Freimuth
!***********************************************************************
CONTAINS
SUBROUTINE wann_anglmom(atoms,usdus,jspin,acof,bcof,ccof, mmn)
USE m_types
IMPLICIT NONE
! .. scalar arguments ..
TYPE(t_atoms),INTENT(in)::atoms
TYPE(t_usdus),INTENT(in)::usdus
INTEGER,INTENT(IN) ::jspin
! .. array arguments ..
COMPLEX, INTENT (in) :: ccof(-atoms%llod:,:,:,:) !ccof(-llod:llod,noccbd,atoms%nlod,natd)
COMPLEX, INTENT (in) :: acof(:,0:,:)!acof(noccbd,0:lmd,natd)
COMPLEX, INTENT (in) :: bcof(:,0:,:)!bcof(noccbd,0:lmd,natd)
COMPLEX, INTENT (inout) :: mmn(:,:,:)!mmn(3,noccbd,noccbd)
! .. local scalars ..
LOGICAL :: l_select
INTEGER :: i,j,l,lo,lop,m,natom,nn,ntyp
INTEGER :: nt1,nt2,lm,n,ll1,indat
COMPLEX :: suma_z,sumb_z
COMPLEX :: suma_p,sumb_p
COMPLEX :: suma_m,sumb_m
COMPLEX :: suma_x,sumb_x
COMPLEX :: suma_y,sumb_y
REAL :: lplus,lminus
! ..
! .. local arrays ..
COMPLEX, ALLOCATABLE :: qlo_z(:,:,:,:,:)
COMPLEX, ALLOCATABLE :: qlo_p(:,:,:,:,:)
COMPLEX, ALLOCATABLE :: qlo_m(:,:,:,:,:)
COMPLEX, ALLOCATABLE :: qaclo_z(:,:,:,:),qbclo_z(:,:,:,:)
COMPLEX, ALLOCATABLE :: qaclo_p(:,:,:,:),qbclo_p(:,:,:,:)
COMPLEX, ALLOCATABLE :: qaclo_m(:,:,:,:),qbclo_m(:,:,:,:)
! ..
! .. intrinsic functions ..
INTRINSIC conjg
ALLOCATE (qlo_z(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%nlod,atoms%ntype) &
,qaclo_z(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%ntype),&
qbclo_z(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%ntype) )
ALLOCATE (qlo_p(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%nlod,atoms%ntype) &
,qaclo_p(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%ntype),&
qbclo_p(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%ntype) )
ALLOCATE (qlo_m(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%nlod,atoms%ntype)&
,qaclo_m(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%ntype),&
qbclo_m(SIZE(acof,1),SIZE(acof,1),atoms%nlod,atoms%ntype) )
INQUIRE(file='select_anglmom',exist=l_select)
WRITE(*,*)'select_anglmom: ',l_select
IF(l_select) THEN
OPEN(866,file='select_anglmom')
READ(866,*)indat
CLOSE(866)
WRITE(*,*)'anglmom for atom=',indat
WRITE(*,*)atoms%ntype
WRITE(*,*)atoms%neq(indat)
ENDIF
!-----> lapw-lapw-Terms
DO i = 1,SIZE(acof,1)
DO j = 1,SIZE(acof,1)
nt1 = 1
DO n = 1,atoms%ntype
nt2 = nt1 + atoms%neq(n) - 1
DO l = 0,atoms%lmax(n)
suma_z = CMPLX(0.,0.); sumb_z = CMPLX(0.,0.)
suma_m = CMPLX(0.,0.); sumb_m = CMPLX(0.,0.)
suma_p = CMPLX(0.,0.); sumb_p = CMPLX(0.,0.)
IF(l_select .AND. (n.NE.indat)) CYCLE
ll1 = l* (l+1)
DO m = -l,l
lm = ll1 + m
lplus=SQRT(REAL( (l-m)*(l+m+1) ) )
lminus=SQRT(REAL( (l+m)*(l-m+1) ) )
DO natom = nt1,nt2
suma_z = suma_z + acof(i,lm,natom)*&
CONJG(acof(j,lm,natom))*REAL(m)
sumb_z = sumb_z + bcof(i,lm,natom)*&
CONJG(bcof(j,lm,natom))*REAL(m)
IF(m+1.LE.l)THEN
suma_p = suma_p + acof(i,lm,natom)*&
CONJG(acof(j,lm+1,natom))*lplus
sumb_p = sumb_p + bcof(i,lm,natom)*&
CONJG(bcof(j,lm+1,natom))*lplus
ENDIF
IF(m-1.GE.-l)THEN
suma_m = suma_m + acof(i,lm,natom)*&
CONJG(acof(j,lm-1,natom))*lminus
sumb_m = sumb_m + bcof(i,lm,natom)*&
CONJG(bcof(j,lm-1,natom))*lminus
ENDIF
ENDDO
ENDDO
mmn(3,j,i) = mmn(3,j,i) + (suma_z+sumb_z*usdus%ddn(l,n,jspin))
suma_x=0.5*(suma_p+suma_m)
sumb_x=0.5*(sumb_p+sumb_m)
mmn(1,j,i) = mmn(1,j,i) + (suma_x+sumb_x*usdus%ddn(l,n,jspin))
suma_y=CMPLX(0.0,-0.5)*(suma_p-suma_m)
sumb_y=CMPLX(0.0,-0.5)*(sumb_p-sumb_m)
mmn(2,j,i) = mmn(2,j,i) + (suma_y+sumb_y*usdus%ddn(l,n,jspin))
ENDDO ! l
nt1 = nt1 + atoms%neq(n)
ENDDO ! n
ENDDO ! j
ENDDO ! i
!---> Terms involving local orbitals.
qlo_z = 0.0; qlo_p = 0.0; qlo_m = 0.0
qaclo_z = 0.0; qaclo_p = 0.0; qaclo_m = 0.0
qbclo_z = 0.0; qbclo_p = 0.0; qbclo_m = 0.0
natom = 0
DO ntyp = 1,atoms%ntype
DO nn = 1,atoms%neq(ntyp)
natom = natom + 1
IF(l_select .AND. (ntyp.NE.indat)) CYCLE
DO lo = 1,atoms%nlo(ntyp)
l = atoms%llo(lo,ntyp)
ll1 = l* (l+1)
DO m = -l,l
lm = ll1 + m
lplus=SQRT(REAL( (l-m)*(l+m+1) ) )
lminus=SQRT(REAL( (l+m)*(l-m+1) ) )
DO i = 1,SIZE(acof,1)
DO j = 1,SIZE(acof,1)
qbclo_z(j,i,lo,ntyp) = qbclo_z(j,i,lo,ntyp) + (&
bcof(i,lm,natom) * CONJG(ccof(m,j,lo