Line data    Source code

       1              :       MODULE m_dylm
       2              :       use m_juDFT
       3              : c.....------------------------------------------------------------------
       4              : c     preparation of dylmt1(=d(ylm)/dtheta),
       5              : c     dylmt2(=d(dylmt1)/dtheta),
       6              : c     dylmf1, dylmf2 are for fai.
       7              : c     dylmtf=d(dylmt1)/df
       8              : c     t.a. june, 1996.
       9              : c.....------------------------------------------------------------------
      10              :       CONTAINS
      11       119768 :       SUBROUTINE dylm3(
      12       119768 :      >                 lmaxd,lmax,v,ylm,
      13       119768 :      <                 dylmt1,dylmt2,dylmf1,dylmf2,dylmtf)
      14              : c
      15              :       use m_constants
      16              :       IMPLICIT NONE
      17              : C     ..
      18              : C     .. Scalar Arguments ..
      19              :       INTEGER, INTENT (IN) :: lmaxd,lmax
      20              : C     ..
      21              : C     .. Array Arguments ..
      22              :       REAL,    INTENT (IN) :: v(3)
      23              :       COMPLEX, INTENT (IN) ::    ylm( (lmaxd+1)**2 )
      24              :       COMPLEX, INTENT (OUT):: dylmt1( (lmaxd+1)**2 )
      25              :       COMPLEX, INTENT (OUT):: dylmt2( (lmaxd+1)**2 )
      26              :       COMPLEX, INTENT (OUT):: dylmf1( (lmaxd+1)**2 )
      27              :       COMPLEX, INTENT (OUT):: dylmf2( (lmaxd+1)**2 )
      28              :       COMPLEX, INTENT (OUT):: dylmtf( (lmaxd+1)**2 )
      29              : C     ..
      30              : C     .. Local Scalars ..
      31              :       INTEGER lm1,lm,lm2,lm1m,lmm1m,lmm,lmm1,lmm2
      32              :       INTEGER l,m,ll1,llm
      33              :       COMPLEX em1f,em2f,ep1f,ep2f
      34              :       REAL cph,rxy,small,sph,x,xy,y
      35              : C     ..
      36              : C     .. Data Statements ..
      37              :       DATA small/1.0e-12/
      38              : c.....------------------------------------------------------------------
      39              : c     ..
      40       119768 :       IF (lmax.GT.lmaxd) THEN
      41            0 :         WRITE (oUnit,*) lmax,lmaxd
      42            0 :          CALL juDFT_error("lmax.GT.lmaxd",calledby="dylm3")
      43              :       ENDIF
      44              : 
      45              : c--->    calculate sin and cos of phi
      46       119768 :       x = v(1)
      47       119768 :       y = v(2)
      48              : 
      49       119768 :       xy = x*x + y*y
      50       119768 :       rxy = sqrt(xy)
      51              : c
      52       119768 :       IF (rxy.gt.small) THEN
      53       119768 :          cph = x/rxy
      54       119768 :          sph = y/rxy
      55              :       ELSE
      56              :          cph = 1.e0
      57              :          sph = 0.e0
      58              :       ENDIF
      59              : 
      60       119768 :       ep1f=cmplx(cph,sph)
      61       119768 :       em1f=conjg(ep1f)
      62       119768 :       ep2f=ep1f*ep1f
      63       119768 :       em2f=em1f*em1f
      64              : c
      65      1289264 :       DO 21 l=0,lmax
      66      1169496 :          ll1 = l*(l+1) 
      67              : 
      68     12826608 :          DO m=-l,l
      69     11657112 :             llm = ll1 + m + 1
      70     11657112 :             dylmt1(llm) = cmplx(0.0,0.0)
      71     12826608 :             dylmt2(llm) = cmplx(0.0,0.0)
      72              :          ENDDO
      73              : 
      74     12826608 :          DO 23 m=-l,l
      75     11657112 :             llm = ll1 + m + 1
      76              : 
      77     11657112 :             lmm1m = l - m - 1
      78     11657112 :             lmm   = l - m
      79     11657112 :             lmm1  = l - m + 1
      80     11657112 :             lmm2  = l - m + 2
      81     11657112 :             lm1m  = l + m - 1
      82     11657112 :             lm    = l + m
      83     11657112 :             lm1   = l + m + 1
      84     11657112 :             lm2   = l + m + 2
      85              : 
      86              :             dylmt2(llm)=dylmt2(llm) -
      87     11657112 :      +                      (lmm*lm1+lmm1*lm)/4.e0*ylm(llm)
      88              : 
      89     11657112 :             IF (m+2.le.l) dylmt2(llm)=dylmt2(llm) +
      90      9437888 :      +           sqrt(real(lmm1m*lmm*lm1*lm2))/4*ylm(llm+2)*em2f
      91              : 
      92     11657112 :             IF (m+1.le.l) dylmt1(llm)=dylmt1(llm) +
      93     10487616 :      +           sqrt(real(lmm*lm1))/2*ylm(llm+1)*em1f
      94              : 
      95     11657112 :             IF (m-1.ge.-l) dylmt1(llm)=dylmt1(llm) -
      96     10487616 :      +           sqrt(real(lm*lmm1))/2*ylm(llm-1)*ep1f
      97              : 
      98     11657112 :             IF (m-2.ge.-l) dylmt2(llm)=dylmt2(llm) +
      99      9437888 :      +           sqrt(real(lmm1*lmm2*lm1m*lm))/4*ylm(llm-2)*ep2f
     100              : 
     101      1169496 :    23    ENDDO
     102              : 
     103     12826608 :          DO m=-l,l
     104     11657112 :             llm = ll1 + m + 1
     105     11657112 :             dylmf1(llm) = ImagUnit * m *    ylm(llm)
     106     11657112 :             dylmf2(llm) = -m * m *    ylm(llm)
     107     12826608 :             dylmtf(llm) = ImagUnit * m * dylmt1(llm)
     108              :          ENDDO
     109              : 
     110       119768 :    21 ENDDO
     111              : 
     112       119768 :       RETURN
     113              :       END SUBROUTINE dylm3
     114              :       END MODULE m_dylm