text:       MODULE m_lhcal
        use m_juDFT
  
  !*********************************************************************
  !     determines the lattice harmonics for the given local
  !     symmetry defined in terms of the rotation matrices.
  !
  !     input:  lmax     max l for lattice harmonics
  !             nrot     number of operations
  !             orth     rotation matrices (cartesian coordinates)
  !                          ( 1,1 1,2 1,3 )( x )   ( x' )
  !                          ( 2,1 2,2 2,3 )( y ) = ( y' )
  !                          ( 3,1 3,2 3,3 )( z )   ( z' )
  !             memd     max. members/lattice harmonic dimensioned for
  !             nlhd     max. number of lattice harmonics dimensioned for
  !
  !     output: nlh      number of harmonic from l=0 to l=lmax
  !             lnu      l value of harmonic
  !             mem      number of members in lattice harmonic
  !             lmnu     lm index of each member of lattice harmonic
  !             c        coefficients of lattice harmonic
  !
  !            m. weinert                           1989/98
  !*********************************************************************
  
        CONTAINS 
        SUBROUTINE lhcal(
       >                 memd,nlhd,lmax,nrot,orth,
       <                 nlh,lnu,mem,lmnu,c)
  !DEC$ NOOPTIMIZE
        USE m_constants
        USE m_gaussp
        USE m_gtest
        USE m_ylm
        IMPLICIT NONE
  
  !---> Arguments
        INTEGER, INTENT (IN) :: memd,nlhd
        INTEGER, INTENT (IN) :: lmax           !--->    max. l to consider
        INTEGER, INTENT (IN) :: nrot           !--->    number of
        REAL,    INTENT (IN) :: orth(3,3,nrot) !--->    rotation matrices
  
        INTEGER, INTENT(OUT) :: nlh
        INTEGER, INTENT(OUT) :: lnu((lmax+1)**2),mem((lmax+1)**2)
        INTEGER, INTENT(OUT) :: lmnu(memd,(lmax+1)**2)
        COMPLEX, INTENT(OUT) ::    c(memd,(lmax+1)**2)
  
  !---> Locals
        INTEGER :: j,l,lh,lm,lm0,lmmax,l2,m,mems,mp,n,nn
        REAL    :: s
        REAL    :: v(3),vr(3),ovlp(0:2*lmax)
        COMPLEX :: a(0:lmax,0:2*lmax,0:lmax)
        COMPLEX, DIMENSION((LMAX+1)**2) :: ylm,ylmr,ylms
  
        INTEGER :: ngpts                                           ! gaussian
        REAL    :: vgauss(3,(2*lmax+1)*(lmax+1 + mod(lmax+1,2)))   ! integration
        REAL    :: wt((2*lmax+1)*( lmax+1 + mod(lmax+1,2) ))       ! points
  
        REAL, PARAMETER :: del = 1.e-5 ! parameter for machine roundoff, etc
  !
  !--->    generate gaussian points and test
  !
  !$OMP MASTER
        ngpts = (2*lmax+1)*( lmax+1 + mod(lmax+1,2) )
        CALL gaussp(
       >            lmax,
       <            vgauss,wt)
        CALL gtest(
       >           lmax,ngpts,vgauss,wt)
  
  !--->    initialize
        lmmax = (lmax+1)**2
        a = cmplx(0.0,0.0)
  !
  !===>    loop over gaussian integration points
  !
        DO nn = 1, ngpts
           v(:) = vgauss(:,nn)
           CALL ylm4(
       >             lmax,v,
       <             ylm)
           ylms(1:lmmax) = ylm(1:lmmax)
  !--->    apply rotations
           DO n =2, nrot
              vr = matmul( orth(:,:,n), v )
              CALL ylm4(
       >                lmax,vr,
       <                ylmr)
              ylms(:) = ylms(:) + ylmr(:)
           ENDDO
           ylms = ylms/nrot
  !--->    obtain coefficients
           DO l = 0, lmax
              lm0 = l*(l+1)+1
              DO mp = 0,l
                 a(mp,0,l) = a(mp,0,l) +
       +                  wt(nn)*conjg(ylm(lm0+mp))*real(ylms(lm0))
              ENDDO
              DO m = 1, l
                 lm = lm0+m
                 DO mp = 0, l
                    a(mp,m,l)   = a(mp,m,l) +
       +                  wt(nn)*conjg(ylm(lm0+mp))* real(ylms(lm))
                    a(mp,m+l,l) = a(mp,m+l,l) +
       +                  wt(nn)*conjg(ylm(lm0+mp))*aimag(ylms(lm))
                 ENDDO
              ENDDO
           ENDDO
        ENDDO
  !
  !===>    orthonormalize the projections for each l (many maybe zero)
  !
        nlh = 0
        DO l = 0, lmax
           lm0 = l*(l+1)+1
           l2  = 2*l
           DO m = 0, l2
  !--->       calculate overlaps with previous harmonics for this l
              DO j = 0, m-1
                 s = real( conjg(a(0,j,l)) * a(0,m,l))
                 DO mp=1,l
                    s = s + 2*real( conjg(a(mp,j,l)) * a(mp,m,l) )
                 ENDDO
                 ovlp(j) = s
              ENDDO
  !--->       Gram-Schmidt
              DO j = 0,m-1
                 a(0:l,m,l) = a(0:l,m,l) - ovlp(j)*a(0:l,j,l)
              ENDDO
  !--->       normalize
              s = real( conjg(a(0,m,l)) * a(0,m,l))
              DO mp = 1, l
                 s = s + 2*real( conjg(a(mp,m,l)) * a(mp,m,l))
              ENDDO
              IF (s.GT.del) THEN
                 s = sqrt(s)
                 a(0:l,m,l) = a(0:l,m,l)/s
  !--->          store lattice harmonic
                 mems = 0
                 nlh = nlh + 1
                 IF (nlh>nlhd+1)  CALL juDFT_error("nlhd too small",
       +              calledby="lhcal")
                 lnu(nlh) = l
                 IF (abs(a(0,m,l)).GT.del) THEN
                    mems = mems+1
                    IF (mems>memd)  CALL juDFT_error("memd too small"
       +                 ,calledby ="lhcal")
                    c(mems,nlh) = a(0,m,l)
                    lmnu(mems,nlh) = lm0
                 ENDIF
                 DO mp=1,l
                    IF( abs(a(mp,m,l)).GT.del) THEN
                       mems = mems + 1
                       IF(mems>memd)  CALL juDFT_error("memd too small"
       +                    ,calledby ="lhcal")
                       c(mems,nlh) = a(mp,m,l)
                       lmnu(mems,nlh) = lm0 + mp
                       mems = mems + 1
                       IF (mems>memd)  CALL juDFT_error("memd too small"
       +                    ,calledby ="lhcal")
                       c(mems,nlh) = ((-1.)**mp)*conjg(c(mems-1,nlh))
                       lmnu(mems,nlh) = lm0 - mp
                    ENDIF
                 ENDDO
                 mem(nlh) = mems
              ELSE
                 a(0:l,m,l) = cmplx(0.0,0.0)
              ENDIF
           ENDDO   ! m = 0, l2
        ENDDO      ! l = 0, lmax
  !
  !===>    test of lattice harmonics using an arbitary point
  !
        v(1) = sqrt(2.0)
        v(2) = sqrt(5.0)
        v(3) = sqrt(17.0)
  !---> generate lattice harmonic for this point
        CALL ylm4(
       >          lmax,v,
       <          ylm)
        DO lh = 1, nlh
           ylms(lh) = cmplx(0.0,0.0)
           DO m = 1, mem(lh)
              ylms(lh) = ylms(lh) + c(m,lh)*ylm(lmnu(m,lh))
           ENDDO
        ENDDO
  !---> rotate point and generate lattice harmonic
        DO n = 2, nrot
           vr = matmul( orth(:,:,n), v )
           CALL ylm4(
       >             lmax,vr,
       <             ylm)
           DO lh = 1, nlh
              ylmr(lh) = cmplx(0.0,0.0)
              DO m = 1, mem(lh)
                 ylmr(lh) = ylmr(lh) + c(m,lh)*ylm(lmnu(m,lh))
              ENDDO
              IF ( abs(ylms(lh)-ylmr(lh)).GT.del ) THEN
                 WRITE (oUnit,'(/," error for operation",i3)') n
                 WRITE (oUnit,'(  " lattice harmonic   ",i3)') lh
                 WRITE (oUnit,'(4f12.6)') ylms(lh),ylmr(lh)
                  CALL juDFT_error("k_lv(Rr)",calledby="lhcal")
              ENDIF
           ENDDO
        ENDDO
  !$OMP END MASTER
        RETURN
        END SUBROUTINE lhcal
        END MODULE m_lhcal
  