Added new wannier routines from JPH

wannier/abcof_small.F

+      MODULE m_abcof_small
+      CONTAINS
+      SUBROUTINE abcof_small(
+     >                 lmaxd,ntypd,neigd,nobd,natd,nop,nvd,jspd,
+     >                 lmd,nbasfcn,llod,nlod,nlotot,invtab,
+     >                 ntype,mrot,ngopr,taual,neq,lmax,rmt,omtil,
+     >                 bmat,bbmat,bkpt,k1,k2,k3,nv,nmat,ne,z,
+     >                 us,dus,uds,duds,ddn,invsat,invsatnr,
+     >                 ulos,uulon,dulon,dulos,llo,nlo,l_dulo,lapw_l,
+     >                 l_noco,l_ss,jspin,alph,beta,qss,kveclo,odi,ods,
+     <                 acof,bcof,ccof,nig)
+c     ************************************************************
+c     subroutine constructs the a,b coefficients of the linearized
+c     m.t. wavefunctions for each band and atom.       c.l. fu
+c     ************************************************************
+#include "cpp_double.h"
+
+      USE m_cotra,     ONLY : cotra3
+      USE m_constants, ONLY : pimach
+      USE m_dotir,     ONLY : dotirp
+      USE m_setabc1locdn
+      USE m_sphbes
+      USE m_dsphbs
+      USE m_abclocdn
+      USE m_ylm
+      USE m_od_types, ONLY : od_inp, od_sym
+
+      IMPLICIT NONE
+C     ..
+C     .. Scalar Arguments ..
+      INTEGER, INTENT (IN) :: lmaxd,ntypd,neigd,nobd,natd,nop,nvd,jspd
+      INTEGER, INTENT (IN) :: ne,ntype,nmat,nbasfcn,llod,nlod,lmd
+      REAL,    INTENT (IN) :: omtil
+      INTEGER, INTENT (IN) :: jspin,nlotot,nig
+      LOGICAL, INTENT (IN) :: l_noco,l_ss
+C     ..
+C     .. Array Arguments ..
+      INTEGER, INTENT (IN) :: mrot(3,3,nop),ngopr(natd),lmax(ntypd)
+      INTEGER, INTENT (IN) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
+      INTEGER, INTENT (IN) :: nv(jspd),lapw_l(ntypd),invtab(nop)
+      INTEGER, INTENT (IN) :: neq(ntypd),invsat(natd),invsatnr(natd)
+      INTEGER, INTENT (IN) :: nlo(ntypd),llo(nlod,ntypd),kveclo(nlotot)
+      REAL,    INTENT (IN) :: bmat(3,3),bbmat(3,3),bkpt(3)
+      REAL,    INTENT (IN) :: taual(3,natd),rmt(ntypd)
+      REAL,    INTENT (IN) :: dus(0:lmaxd,ntypd),duds(0:lmaxd,ntypd)
+      REAL,    INTENT (IN) ::  us(0:lmaxd,ntypd), uds(0:lmaxd,ntypd)
+      REAL,    INTENT (IN) :: ulos(nlod,ntypd),uulon(nlod,ntypd)
+      REAL,    INTENT (IN) :: dulon(nlod,ntypd),dulos(nlod,ntypd)
+      REAL,    INTENT (IN) :: ddn(0:lmaxd,ntypd)
+      REAL,    INTENT (IN) :: alph(ntypd),beta(ntypd),qss(3)
+#if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )
+      COMPLEX, INTENT (IN) :: z(nbasfcn,neigd)
+#else
+      REAL,    INTENT (IN) :: z(nbasfcn,neigd)
+#endif
+      LOGICAL, INTENT (IN) :: l_dulo(nlod,ntypd)
+      COMPLEX, INTENT (OUT):: acof(nobd,0:lmd,natd)
+      COMPLEX, INTENT (OUT):: bcof(nobd,0:lmd,natd)
+      COMPLEX, INTENT (OUT):: ccof(-llod:llod,nobd,nlod,natd)
+c-odim
+      TYPE (od_inp), INTENT (IN) :: odi
+      TYPE (od_sym), INTENT (IN) :: ods
+c+odim
+C     ..
+C     .. Local Scalars ..
+      COMPLEX phase,cexp,c_0,c_1,c_2,ci
+      REAL const,df,r1,s,tmk,wronk,tpi,qss1,qss2,qss3
+      INTEGER i,j,k,l,ll1,lm,m,n,nap,natom,nn,iatom,jatom,lmp
+      INTEGER inv_f,ie,ilo,kspin,iintsp,nintsp,nvmax,lo,inap
+C     ..
+C     .. Local Arrays ..
+      INTEGER kvec(2*(2*llod+1),nlod,natd)
+      INTEGER nbasf0(nlod,natd),nkvec(nlod,natd)
+      REAL dfj(0:lmaxd),fj(0:lmaxd),fk(3),fkp(3),fkr(3)
+      REAL alo1(nlod,ntypd),blo1(nlod,ntypd),clo1(nlod,ntypd)
+      COMPLEX ylm( (lmaxd+1)**2 )
+      COMPLEX ccchi(2,2)
+      LOGICAL enough(natd),apw(0:lmaxd,ntypd)
+#if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )
+      COMPLEX, ALLOCATABLE :: work(:)
+#else
+      REAL,    ALLOCATABLE :: work(:)
+#endif
+C     ..
+C     .. Intrinsic Functions ..
+      INTRINSIC cmplx,conjg,exp,sqrt
+C     ..
+      ci = cmplx(0.0,1.0)
+      tpi = 2 * pimach()
+      const = 2 * tpi/sqrt(omtil)
+c
+      acof(:,:,:) = cmplx(0.0,0.0)
+      bcof(:,:,:) = cmplx(0.0,0.0)
+      ALLOCATE ( work(nobd) )
+C     ..
+c+APW_LO
+      DO n = 1, ntype
+         DO l = 0,lmax(n)
+           apw(l,n) = .false.
+           DO lo = 1,nlo(n)
+             IF (l_dulo(lo,n)) apw(l,n) = .true.
+           ENDDO
+#ifdef CPP_APW
+           IF (lapw_l(n).GE.l) apw(l,n) = .false.
+#endif
+         ENDDO
+         DO lo = 1,nlo(n)
+           IF (l_dulo(lo,n)) apw(llo(lo,n),n) = .true.
+         ENDDO
+      ENDDO
+c+APW_LO
+c
+      nintsp = 1
+      IF (l_ss) nintsp = 2
+c---> loop over the interstitial spin
+      DO iintsp = 1,nintsp
+      nvmax = nv(jspin)
+      IF (l_ss) nvmax = nv(iintsp)
+c
+      CALL setabc1locdn(
+     >                  ntypd,natd,nlod,llod,nobd,lmaxd,nvmax,
+     >                  nmat,ne,ntype,neq,l_noco,iintsp,
+     >                  nlo,llo,l_dulo,invsat,invsatnr,ddn,
+     >                  us,dus,uds,duds,dulos,ulos,dulon,uulon,
+     >                  nlotot,kveclo,
+     <                  enough,nkvec,kvec,nbasf0,ccof,
+     <                  alo1,blo1,clo1)
+c
+      IF (iintsp .EQ. 1) THEN
+         qss1= - qss(1)/2
+         qss2= - qss(2)/2
+         qss3= - qss(3)/2
+      ELSE
+         qss1= + qss(1)/2
+         qss2= + qss(2)/2
+         qss3= + qss(3)/2
+      ENDIF
+
+      k = nig
+
+c---> loop over atom types
+      natom = 0
+      DO n = 1,ntype
+c  ----> loop over equivalent atoms
+         DO nn = 1,neq(n)
+            natom = natom + 1
+           IF ((invsat(natom).EQ.0) .OR. (invsat(natom).EQ.1)) THEN
+c--->        loop over lapws
+             !DO k = 1,nvmax
+               IF (.NOT.l_noco) THEN
+                 DO i = 1,ne
+                   work(i) = z(k,i)
+                 ENDDO
+               ENDIF
+
+               IF (l_noco) THEN
+c--->            generate the complex conjgates of the spinors (chi)
+                 ccchi(1,1) = conjg( exp(-ci*alph(n)/2)*cos(beta(n)/2))
+                 ccchi(1,2) = conjg(-exp(-ci*alph(n)/2)*sin(beta(n)/2))
+                 ccchi(2,1) = conjg( exp( ci*alph(n)/2)*sin(beta(n)/2))
+                 ccchi(2,2) = conjg( exp( ci*alph(n)/2)*cos(beta(n)/2))
+                 IF (l_ss) THEN
+c--->              the coefficients of the spin-down basis functions are
+c--->              stored in the second half of the eigenvector
+                   kspin = (iintsp-1)*(nv(1)+nlotot)
+                   DO i = 1,ne
+                      work(i) = ccchi(iintsp,jspin)*z(kspin+k,i)
+                   ENDDO
+                 ELSE
+c--->              perform sum over the two interstitial spin directions
+c--->              and take into account the spin boundary conditions
+c--->              (jspin counts the local spin directions inside each MT)
+                   kspin = nv(1)+nlotot
+                   DO i = 1,ne
+                      work(i) = ccchi(1,jspin)*z(k,i)
+     +                        + ccchi(2,jspin)*z(kspin+k,i)
+                   ENDDO
+                 ENDIF
+               ENDIF ! (l_noco)
+               IF (l_ss) THEN
+                 fk(1) = bkpt(1) + k1(k,iintsp) + qss1
+                 fk(2) = bkpt(2) + k2(k,iintsp) + qss2
+                 fk(3) = bkpt(3) + k3(k,iintsp) + qss3
+               ELSE
+                 fk(1) = bkpt(1) + k1(k,jspin) + qss1
+                 fk(2) = bkpt(2) + k2(k,jspin) + qss2
+                 fk(3) = bkpt(3) + k3(k,jspin) + qss3
+               ENDIF ! (l_ss)
+               s = dotirp(fk,fk,bbmat)
+               s = sqrt(s)
+               r1 = rmt(n)*s
+               CALL sphbes(
+     >                     lmax(n),r1,
+     <                     fj)
+               CALL dsphbs(
+     >                     lmax(n),r1,fj,
+     <                     dfj)
+c  ----> construct a and b coefficients
+               DO l = 0,lmax(n)
+                 df = s*dfj(l)
+                 wronk = uds(l,n)*dus(l,n) - us(l,n)*duds(l,n)
+                 IF (apw(l,n)) THEN
+                   fj(l) = 1.0*const * fj(l)/us(l,n)
+                   dfj(l) = 0.0d0
+                 ELSE
+                   dfj(l) = const* (dus(l,n)*fj(l)-df*us(l,n))/wronk
+                   fj(l) = const* (df*uds(l,n)-fj(l)*duds(l,n))/wronk
+                 ENDIF
+               ENDDO ! loop over l
+               tmk = tpi* (fk(1)*taual(1,natom)+
+     +                     fk(2)*taual(2,natom)+
+     +                     fk(3)*taual(3,natom))
+               phase = cmplx(cos(tmk),sin(tmk))
+               IF (odi%d1) THEN
+                  inap = ods%ngopr(natom)
+               ELSE
+                  nap = ngopr(natom)
+                  inap = invtab(nap)
+               END IF
+               DO j = 1,3
+                 fkr(j) = 0.
+                 DO i = 1,3
+                   IF (odi%d1) THEN
+                     fkr(j) = fkr(j) + fk(i)*ods%mrot(i,j,inap)
+                   ELSE
+                     fkr(j) = fkr(j) + fk(i)*mrot(i,j,inap)
+                   END IF
+                 ENDDO
+               ENDDO
+               CALL cotra3(fkr,fkp,bmat)
+c     ----> generate spherical harmonics
+               CALL ylm4(
+     >                   lmax(n),fkp,
+     <                   ylm)
+c     ----> loop over l
+               DO l = 0,lmax(n)
+                 ll1 = l* (l+1)
+c     ----> loop over m
+                 DO m = -l,l
+                   lm = ll1 + m
+                   c_0 = conjg(ylm(lm+1))*phase
+                   c_1 = c_0 *  fj(l)
+                   c_2 = c_0 * dfj(l)
+c     ----> loop over bands
+                   DO i = 1,ne
+                     acof(i,lm,natom) = acof(i,lm,natom) + 
+     +                                  c_1 * work(i)
+                   ENDDO
+                   DO i = 1,ne
+                     bcof(i,lm,natom) = bcof(i,lm,natom) +
+     +                                  c_2 * work(i)
+                   ENDDO
+#if ( defined(CPP_SOC) && defined(CPP_INVERSION) )
+                   IF (invsat(natom).EQ.1) THEN
+                     jatom = invsatnr(natom)
+                     lmp = ll1 - m
+                     inv_f = (-1)**(l-m)
+                     c_1 =  conjg(c_1) * inv_f
+                     c_2 =  conjg(c_2) * inv_f
+                     CALL CPP_BLAS_caxpy(ne,c_1,work,1,
+     X                                   acof(1,lmp,jatom),1)
+                     CALL CPP_BLAS_caxpy(ne,c_2,work,1,
+     X                                   bcof(1,lmp,jatom),1)
+                   ENDIF
+#endif
+                 ENDDO ! loop over m
+               ENDDO ! loop over l
+               IF (.NOT.enough(natom)) THEN
+                 write(*,*)'.not.enough(natom)'
+                 CALL abclocdn(
+     >                 nobd,natd,nlod,llod,lmaxd,neigd,ntypd,
+     >                 lmd,nbasfcn,nlo,llo,invsat,invsatnr,
+     >                 l_noco,ccchi(1,jspin),kspin,l_ss,iintsp,
+     >                 const,rmt(n),phase,ylm,n,natom,k,s,nvmax,
+     >                 ne,z,nbasf0,alo1,blo1,clo1,kvec(1,1,natom),
+     <                 nkvec,enough,acof,bcof,ccof)
+               ENDIF
+             !ENDDO ! loop over LAPWs
+           ENDIF  ! invsatom == ( 0 v 1 )
+         ENDDO    ! loop over equivalent atoms
+      ENDDO       ! loop over atom types
+      ENDDO       ! loop over interstitial spin
+
+#if ( defined(CPP_SOC) && defined(CPP_INVERSION) )
+!
+!                           -p,n       (l+m)   p,n  *
+! Usually, we exploit that A     = (-1)      (A    )  if p and -p are the positions
+!                           l,m                l,-m  
+! of two atoms related by inversion symmetry and the coefficients are considered to
+! be in the local frame of the representative atom. This is possible, if z is real.
+! After SOC, however, the eigenvectors z are complex and this is no longer possible
+! so the z has to enter, not z*. This is done within the k-loop.
+!                                    -p,n       m   p,n  *
+! When called from hsohelp, we need A     = (-1)  (A    ) because we don't have to
+!                                     l,m           l,-m                    l 
+! rotate, but in the sums in hsoham only products A*  A   enter and the (-1) cancels.
+!                                                  lm  lm
+#else
+      iatom = 0
+      DO n = 1,ntype
+         DO nn = 1,neq(n)
+            iatom = iatom + 1
+            IF (invsat(iatom).EQ.1) THEN
+              jatom = invsatnr(iatom)
+              cexp = exp(tpi*ci*dot_product(taual(:,jatom)
+     +             + taual(:,iatom),(/bkpt(1),bkpt(2),bkpt(3)/)))
+               DO ilo = 1,nlo(n)
+                  l = llo(ilo,n)
+                  DO m = -l,l
+                     inv_f = (-1.0)**(m+l)
+                     DO ie = 1,ne
+                        ccof(m,ie,ilo,jatom) = inv_f * cexp *
+     +                               conjg(  ccof(-m,ie,ilo,iatom))
+                     ENDDO
+                  ENDDO
+               ENDDO
+               DO l = 0,lmax(n)
+                  ll1 = l* (l+1)
+                  DO m =-l,l
+                     lm  = ll1 + m
+                     lmp = ll1 - m
+                     inv_f = (-1.0)**(m+l)
+                     DO ie = 1,ne
+                        acof(ie,lm,jatom) = inv_f * cexp * 
+     *                                      conjg(acof(ie,lmp,iatom))
+                     ENDDO
+                     DO ie = 1,ne
+                        bcof(ie,lm,jatom) = inv_f * cexp * 
+     *                                      conjg(bcof(ie,lmp,iatom))
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDIF
+         ENDDO
+      ENDDO
+#endif
+      DEALLOCATE ( work )
+c
+      END SUBROUTINE abcof_small
+      END MODULE m_abcof_small

wannier/uhu/d2fdz2.F

+      MODULE m_d2fdz2
+      CONTAINS
+
+      SUBROUTINE d2fdz2(jmtd,jri,rmsh,dx,f,fac,d2f)
+
+      IMPLICIT NONE
+
+      REAL difcub
+      EXTERNAL difcub
+      
+      REAL, INTENT(INOUT) :: d2f(jmtd)
+      REAL, INTENT(IN)    :: fac(jmtd)
+      REAL, INTENT(IN)    :: f(jmtd)
+      REAL, INTENT(IN)    :: rmsh(jmtd)
+      REAL, INTENT(IN)    :: dx
+
+      INTEGER, INTENT(IN) :: jri
+      INTEGER, INTENT(IN) :: jmtd
+
+      REAL, ALLOCATABLE :: fr(:),dfr(:)
+      INTEGER :: i
+
+      allocate( dfr(jri),fr(jri) )
+      DO i=1,jri
+       fr(i) = f(i)*fac(i)
+      ENDDO
+
+      dfr(1) = difcub( rmsh(1),fr(1),rmsh(1) )
+      DO i = 2, jri-2
+         dfr(i) = difcub( rmsh(i-1),fr(i-1),rmsh(i) )
+      ENDDO
+      dfr(jri-1) = difcub( rmsh(jri-3),fr(jri-3),rmsh(jri-1) )
+      dfr(jri) = difcub( rmsh(jri-3),fr(jri-3),rmsh(jri) )
+
+
+      d2f(1) = difcub( rmsh(1),dfr(1),rmsh(1) )
+      DO i = 2, jri-2
+         d2f(i) = difcub( rmsh(i-1),dfr(i-1),rmsh(i) )
+      ENDDO
+      d2f(jri-1) = difcub( rmsh(jri-3),dfr(jri-3),rmsh(jri-1) )
+      d2f(jri) = difcub( rmsh(jri-3),dfr(jri-3),rmsh(jri) )
+
+
+      deallocate( dfr,fr )
+
+      END SUBROUTINE d2fdz2
+      END MODULE m_d2fdz2

wannier/uhu/d2fdz2cmplx.F

+      MODULE m_d2fdz2cmplx
+      CONTAINS
+
+      SUBROUTINE d2fdz2cmplx(jmtd,jri,rmsh,dx,f,fac,d2f)
+
+      IMPLICIT NONE
+
+      REAL difcub
+      EXTERNAL difcub
+      
+      COMPLEX, INTENT(INOUT) :: d2f(jmtd)
+      COMPLEX, INTENT(IN)    :: fac(jmtd)
+
+      REAL, INTENT(IN)    :: f(jmtd)
+      REAL, INTENT(IN)    :: rmsh(jmtd)
+      REAL, INTENT(IN)    :: dx
+
+      INTEGER, INTENT(IN) :: jri
+      INTEGER, INTENT(IN) :: jmtd
+
+      REAL, ALLOCATABLE :: fr(:),fi(:),dfr(:),dfi(:)
+      INTEGER :: i
+
+      allocate( dfr(jri),dfi(jri),fr(jri),fi(jri) )
+      DO i=1,jri
+       fr(i) = f(i)*real(fac(i))
+       fi(i) = f(i)*aimag(fac(i))
+      ENDDO
+
+      dfr(1) = difcub( rmsh(1),fr(1),rmsh(1) )
+      dfi(1) = difcub( rmsh(1),fi(1),rmsh(1) )
+      DO i = 2, jri-2
+         dfr(i) = difcub( rmsh(i-1),fr(i-1),rmsh(i) )
+         dfi(i) = difcub( rmsh(i-1),fi(i-1),rmsh(i) )
+      ENDDO
+      dfr(jri-1) = difcub( rmsh(jri-3),fr(jri-3),rmsh(jri-1) )
+      dfi(jri-1) = difcub( rmsh(jri-3),fi(jri-3),rmsh(jri-1) )
+      dfr(jri) = difcub( rmsh(jri-3),fr(jri-3),rmsh(jri) )
+      dfi(jri) = difcub( rmsh(jri-3),fi(jri-3),rmsh(jri) )
+
+
+      d2f(1) = cmplx( difcub( rmsh(1),dfr(1),rmsh(1) ),
+     >                difcub( rmsh(1),dfi(1),rmsh(1) ) )
+      DO i = 2, jri-2
+         d2f(i) = cmplx( difcub( rmsh(i-1),dfr(i-1),rmsh(i) ),
+     >                   difcub( rmsh(i-1),dfi(i-1),rmsh(i) ) )
+      ENDDO
+      d2f(jri-1) = cmplx( difcub( rmsh(jri-3),dfr(jri-3),rmsh(jri-1) ),
+     >                    difcub( rmsh(jri-3),dfi(jri-3),rmsh(jri-1) ) )
+      d2f(jri) = cmplx( difcub( rmsh(jri-3),dfr(jri-3),rmsh(jri) ),
+     >                  difcub( rmsh(jri-3),dfi(jri-3),rmsh(jri) ) )
+
+      deallocate( dfr,dfi,fr,fi )
+
+c      d2f = cmplx(0.,0.)
+c      d2f(1) = (f(3)-2*f(2)+f(1))/dx/dx
+c      do i=2,jri-1
+c       d2f(i) = (f(i+1)-2*f(i)+f(i-1))/dx/dx
+c      enddo
+c      d2f(jri) = (f(jri-2)-2*f(jri-1)+f(jri))/dx/dx 
+
+      END SUBROUTINE d2fdz2cmplx
+      END MODULE m_d2fdz2cmplx

wannier/uhu/dujdr.F

+      MODULE m_dujdr
+      CONTAINS
+
+      SUBROUTINE dujdr(jmtd,jri,rmsh,dx,ub,j2,b,l,lmaxd,dub)
+
+      USE m_constants
+      IMPLICIT NONE
+
+      REAL difcub
+      EXTERNAL difcub
+      
+      REAL, INTENT(IN) :: rmsh(jmtd),dx
+      REAL, INTENT(IN) :: ub(jmtd,2)       ! u(b2)
+      REAL, INTENT(IN) :: j2(0:lmaxd,jmtd) ! j_l(b2*r)
+      REAL, INTENT(IN) :: b                ! b2
+      INTEGER, INTENT(IN) :: jri,jmtd,lmaxd
+      INTEGER, INTENT(IN) :: l           ! l of sph. Bessel j2
+
+      REAL, INTENT(OUT) :: dub(jmtd,2)
+      REAL :: xi,t(jri,2)
+      INTEGER :: i,j
+
+      ! derivatives d/dr for large and small component of q
+      DO i=1,jri
+       t(i,:) = ub(i,:) / rmsh(i) * j2(l,i)
+      ENDDO
+
+      DO j = 1, 2
+         ! derivative at 1st point
+         dub(1,j) = difcub( rmsh(1),t(1,j),rmsh(1) )
+
+         ! derivative at 2nd...(jri-2)th point
+         DO i = 2, jri-2
+            dub(i,j) = difcub( rmsh(i-1),t(i-1,j),rmsh(i) )
+         ENDDO
+
+         ! derivative at last two points
+         dub(jri-1,j) = difcub( rmsh(jri-3),t(jri-3,j),rmsh(jri-1) )
+         dub(jri,j) = difcub( rmsh(jri-3),t(jri-3,j),rmsh(jri) )
+      ENDDO
+
+      DO i=1,jri
+       dub(i,:) = dub(i,:)*rmsh(i)
+      ENDDO
+
+      ! complete d/dr (ub*j2) = ub'j2 + ub j2' with sph. Bessel func. j
+      ! rule:  j'_{l}(ar) = a*j_{l-1}(ar) - (l+1)/r*j_{l}(ar)
+c      IF(l.ne.0) THEN
+c       DO i=1,jri
+c        xi = rmsh(i)
+c        dub(i,:) = dub(i,:) * j2(l,i) * xi
+c     >           + ub(i,:) *( j2(l-1,i)*b - (l+1)/xi*j2(l,i) ) 
+c       ENDDO
+c      ELSE
+c       DO i=1,jri
+c        xi = rmsh(i)
+c        dub(i,:) = dub(i,:) * j2(l,i) * xi
+c     >           - ub(i,:) * j2(1,i) * b
+c       ENDDO
+c      ENDIF
+      
+      END SUBROUTINE dujdr
+      END MODULE m_dujdr

wannier/uhu/wann_uHu.F

+c*******************************************c
+c   Set up uHu matrix necessary for         c
+c   Wannier-based calc. of orbital moment   c
+c*******************************************c
+c   keyword is 'matrixuhu' in wann_inp      c
+c*******************************************c
+c   uHu = < u_{k+b1} | H_{k} | u_{k+b2} >   c
+c                                           c
+c   Contributions to Hamiltonian:           c
+c       (i)   interstitial                  c
+c       (ii)  muffin tin  (a) spherical     c
+c                         (b) non-sph.      c
+c                         (c) SOC           c
+c       (iii) vacuum                        c
+c*******************************************c
+c                  J.-P. Hanke, Dec. 2015   c
+c*******************************************c
+      MODULE m_wann_uHu
+      USE m_fleurenv
+      CONTAINS
+      SUBROUTINE wann_uHu(
+     >      l_dulo,l_noco,l_ss,lmaxd,ntypd,
+     >      neigd,natd,nop,nvd,jspd,nbasfcn,llod,nlod,ntype,
+     >      nwdd,omtil,nlo,llo,lapw_l,invtab,mrot,ngopr,neq,lmax,
+     >      invsat,invsatnr,nkpt,taual,rmt,amat,bmat,bbmat,alph,
+     >      beta,qss,sk2,phi2,odi,ods,irank,isize,n3d,nmzxyd,nmzd,
+     >      jmtd,nlhd,nq3,nvac,invs,invs2,film,nlh,jri,ntypsd,
+     >      ntypsy,jspins,nkptd,dx,n2d,rmsh,e1s,e2s,ulo_der,
+     >      ustep,ig,k1d,k2d,k3d,rgphs,slice,kk,nnne,
+     >      z1,nv2d,nmzxy,nmz,delz,zrfs,ig2,area,tau,zatom,nq2,kv2,nop2,
+     >      volint,symor,pos,ef,irecl,l_soc,
+     >      memd,lnonsph,clnu,lmplmd,mlh,nmem,llh,lo1l,
+     >      theta,phi,soc_opt,
+     >      l_ms,l_sgwf,l_socgwf,aux_latt_const,
+     >      param_file,param_vec,nparampts,param_alpha,l_dim)
+
+      use m_types
+      use m_wann_mmnk_symm
+      use m_wann_rw_eig
+      use m_abcof
+      use m_radfun
+      use m_radflo
+      use m_cdnread, only : cdn_read0, cdn_read
+      use m_od_types, only : od_inp, od_sym
+      use m_loddop
+      use m_constants, only : pimach
+      use m_wann_projmethod
+      use m_wann_abinv
+      use m_wann_kptsrotate
+      use m_wann_read_inp
+      use m_matmul,only : matmul3,matmul3r