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      122680 :       SUBROUTINE dylm3(
      12      122680 :      >                 lmaxd,lmax,v,ylm,
      13      122680 :      <                 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      122680 :       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      122680 :       x = v(1)
      47      122680 :       y = v(2)
      48             : 
      49      122680 :       xy = x*x + y*y
      50      122680 :       rxy = sqrt(xy)
      51             : c
      52      122680 :       IF (rxy.gt.small) THEN
      53      122680 :          cph = x/rxy
      54      122680 :          sph = y/rxy
      55             :       ELSE
      56             :          cph = 1.e0
      57             :          sph = 0.e0
      58             :       ENDIF
      59             : 
      60      122680 :       ep1f=cmplx(cph,sph)
      61      122680 :       em1f=conjg(ep1f)
      62      122680 :       ep2f=ep1f*ep1f
      63      122680 :       em2f=em1f*em1f
      64             : c
      65     1312560 :       DO 21 l=0,lmax
      66     1189880 :          ll1 = l*(l+1) 
      67             : 
      68    12989680 :          DO m=-l,l
      69    11799800 :             llm = ll1 + m + 1
      70    11799800 :             dylmt1(llm) = cmplx(0.0,0.0)
      71    12989680 :             dylmt2(llm) = cmplx(0.0,0.0)
      72             :          ENDDO
      73             : 
      74    12989680 :          DO 23 m=-l,l
      75    11799800 :             llm = ll1 + m + 1
      76             : 
      77    11799800 :             lmm1m = l - m - 1
      78    11799800 :             lmm   = l - m
      79    11799800 :             lmm1  = l - m + 1
      80    11799800 :             lmm2  = l - m + 2
      81    11799800 :             lm1m  = l + m - 1
      82    11799800 :             lm    = l + m
      83    11799800 :             lm1   = l + m + 1
      84    11799800 :             lm2   = l + m + 2
      85             : 
      86             :             dylmt2(llm)=dylmt2(llm) -
      87    11799800 :      +                      (lmm*lm1+lmm1*lm)/4.e0*ylm(llm)
      88             : 
      89    11799800 :             IF (m+2.le.l) dylmt2(llm)=dylmt2(llm) +
      90     9542720 :      +           sqrt(real(lmm1m*lmm*lm1*lm2))/4*ylm(llm+2)*em2f
      91             : 
      92    11799800 :             IF (m+1.le.l) dylmt1(llm)=dylmt1(llm) +
      93    10609920 :      +           sqrt(real(lmm*lm1))/2*ylm(llm+1)*em1f
      94             : 
      95    11799800 :             IF (m-1.ge.-l) dylmt1(llm)=dylmt1(llm) -
      96    10609920 :      +           sqrt(real(lm*lmm1))/2*ylm(llm-1)*ep1f
      97             : 
      98    11799800 :             IF (m-2.ge.-l) dylmt2(llm)=dylmt2(llm) +
      99     9542720 :      +           sqrt(real(lmm1*lmm2*lm1m*lm))/4*ylm(llm-2)*ep2f
     100             : 
     101     1189880 :    23    ENDDO
     102             : 
     103    12989680 :          DO m=-l,l
     104    11799800 :             llm = ll1 + m + 1
     105    11799800 :             dylmf1(llm) = ImagUnit * m *    ylm(llm)
     106    11799800 :             dylmf2(llm) = -m * m *    ylm(llm)
     107    12989680 :             dylmtf(llm) = ImagUnit * m * dylmt1(llm)
     108             :          ENDDO
     109             : 
     110      122680 :    21 ENDDO
     111             : 
     112      122680 :       RETURN
     113             :       END SUBROUTINE dylm3
     114             :       END MODULE m_dylm