LCOV - code coverage report
Current view: top level - init - lhcal.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 88 96 91.7 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_lhcal
       2             :       use m_juDFT
       3             : 
       4             : !*********************************************************************
       5             : !     determines the lattice harmonics for the given local
       6             : !     symmetry defined in terms of the rotation matrices.
       7             : !
       8             : !     input:  lmax     max l for lattice harmonics
       9             : !             nrot     number of operations
      10             : !             orth     rotation matrices (cartesian coordinates)
      11             : !                          ( 1,1 1,2 1,3 )( x )   ( x' )
      12             : !                          ( 2,1 2,2 2,3 )( y ) = ( y' )
      13             : !                          ( 3,1 3,2 3,3 )( z )   ( z' )
      14             : !             memd     max. members/lattice harmonic dimensioned for
      15             : !             nlhd     max. number of lattice harmonics dimensioned for
      16             : !
      17             : !     output: nlh      number of harmonic from l=0 to l=lmax
      18             : !             lnu      l value of harmonic
      19             : !             mem      number of members in lattice harmonic
      20             : !             lmnu     lm index of each member of lattice harmonic
      21             : !             c        coefficients of lattice harmonic
      22             : !
      23             : !            m. weinert                           1989/98
      24             : !*********************************************************************
      25             : 
      26             :       CONTAINS 
      27         108 :       SUBROUTINE lhcal(
      28         108 :      >                 memd,nlhd,lmax,nrot,orth,
      29         108 :      <                 nlh,lnu,mem,lmnu,c)
      30             : !DEC$ NOOPTIMIZE
      31             :       USE m_gaussp
      32             :       USE m_gtest
      33             :       USE m_ylm
      34             :       IMPLICIT NONE
      35             : 
      36             : !---> Arguments
      37             :       INTEGER, INTENT (IN) :: memd,nlhd
      38             :       INTEGER, INTENT (IN) :: lmax           !--->    max. l to consider
      39             :       INTEGER, INTENT (IN) :: nrot           !--->    number of
      40             :       REAL,    INTENT (IN) :: orth(3,3,nrot) !--->    rotation matrices
      41             : 
      42             :       INTEGER, INTENT(OUT) :: nlh
      43             :       INTEGER, INTENT(OUT) :: lnu((lmax+1)**2),mem((lmax+1)**2)
      44             :       INTEGER, INTENT(OUT) :: lmnu(memd,(lmax+1)**2)
      45             :       COMPLEX, INTENT(OUT) ::    c(memd,(lmax+1)**2)
      46             : 
      47             : !---> Locals
      48             :       INTEGER :: j,l,lh,lm,lm0,lmmax,l2,m,mems,mp,n,nn
      49             :       REAL    :: s
      50         216 :       REAL    :: v(3),vr(3),ovlp(0:2*lmax)
      51         216 :       COMPLEX :: a(0:lmax,0:2*lmax,0:lmax)
      52         216 :       COMPLEX, DIMENSION((LMAX+1)**2) :: ylm,ylmr,ylms
      53             : 
      54             :       INTEGER :: ngpts                                           ! gaussian
      55         216 :       REAL    :: vgauss(3,(2*lmax+1)*(lmax+1 + mod(lmax+1,2)))   ! integration
      56         216 :       REAL    :: wt((2*lmax+1)*( lmax+1 + mod(lmax+1,2) ))       ! points
      57             : 
      58             :       REAL, PARAMETER :: del = 1.e-5 ! parameter for machine roundoff, etc
      59             : !
      60             : !--->    generate gaussian points and test
      61             : !
      62         108 : !$OMP MASTER
      63         108 :       ngpts = (2*lmax+1)*( lmax+1 + mod(lmax+1,2) )
      64             :       CALL gaussp(
      65             :      >            lmax,
      66         108 :      <            vgauss,wt)
      67             :       CALL gtest(
      68         108 :      >           lmax,ngpts,vgauss,wt)
      69             : 
      70             : !--->    initialize
      71         108 :       lmmax = (lmax+1)**2
      72        1172 :       a = cmplx(0.0,0.0)
      73             : !
      74             : !===>    loop over gaussian integration points
      75             : !
      76       45972 :       DO nn = 1, ngpts
      77       22932 :          v(:) = vgauss(:,nn)
      78             :          CALL ylm4(
      79             :      >             lmax,v,
      80       22932 :      <             ylm)
      81     2798344 :          ylms(1:lmmax) = ylm(1:lmmax)
      82             : !--->    apply rotations
      83      242544 :          DO n =2, nrot
      84      219612 :             vr = matmul( orth(:,:,n), v )
      85             :             CALL ylm4(
      86             :      >                lmax,vr,
      87      219612 :      <                ylmr)
      88    29154796 :             ylms(:) = ylms(:) + ylmr(:)
      89             :          ENDDO
      90       22932 :          ylms = ylms/nrot
      91             : !--->    obtain coefficients
      92      270292 :          DO l = 0, lmax
      93      247252 :             lm0 = l*(l+1)+1
      94     1758584 :             DO mp = 0,l
      95             :                a(mp,0,l) = a(mp,0,l) +
      96     1758584 :      +                  wt(nn)*conjg(ylm(lm0+mp))*real(ylms(lm0))
      97             :             ENDDO
      98     2775412 :             DO m = 1, l
      99     1264080 :                lm = lm0+m
     100    12175172 :                DO mp = 0, l
     101             :                   a(mp,m,l)   = a(mp,m,l) +
     102    10663840 :      +                  wt(nn)*conjg(ylm(lm0+mp))* real(ylms(lm))
     103             :                   a(mp,m+l,l) = a(mp,m+l,l) +
     104    11927920 :      +                  wt(nn)*conjg(ylm(lm0+mp))*aimag(ylms(lm))
     105             :                ENDDO
     106             :             ENDDO
     107             :          ENDDO
     108             :       ENDDO
     109             : !
     110             : !===>    orthonormalize the projections for each l (many maybe zero)
     111             : !
     112         108 :       nlh = 0
     113         108 :       DO l = 0, lmax
     114        1064 :          lm0 = l*(l+1)+1
     115        1064 :          l2  = 2*l
     116       12052 :          DO m = 0, l2
     117             : !--->       calculate overlaps with previous harmonics for this l
     118      157864 :             DO j = 0, m-1
     119       73438 :                s = real( conjg(a(0,j,l)) * a(0,m,l))
     120      661440 :                DO mp=1,l
     121      661440 :                   s = s + 2*real( conjg(a(mp,j,l)) * a(mp,m,l) )
     122             :                ENDDO
     123       84426 :                ovlp(j) = s
     124             :             ENDDO
     125             : !--->       Gram-Schmidt
     126      157864 :             DO j = 0,m-1
     127       84426 :                a(0:l,m,l) = a(0:l,m,l) - ovlp(j)*a(0:l,j,l)
     128             :             ENDDO
     129             : !--->       normalize
     130       10988 :             s = real( conjg(a(0,m,l)) * a(0,m,l))
     131       84426 :             DO mp = 1, l
     132       84426 :                s = s + 2*real( conjg(a(mp,m,l)) * a(mp,m,l))
     133             :             ENDDO
     134       12052 :             IF (s.GT.del) THEN
     135        2228 :                s = sqrt(s)
     136        2228 :                a(0:l,m,l) = a(0:l,m,l)/s
     137             : !--->          store lattice harmonic
     138        2228 :                mems = 0
     139        2228 :                nlh = nlh + 1
     140        2228 :                IF (nlh>nlhd+1)  CALL juDFT_error("nlhd too small",
     141           0 :      +              calledby="lhcal")
     142        2228 :                lnu(nlh) = l
     143        2228 :                IF (abs(a(0,m,l)).GT.del) THEN
     144         818 :                   mems = mems+1
     145         818 :                   IF (mems>memd)  CALL juDFT_error("memd too small"
     146           0 :      +                 ,calledby ="lhcal")
     147         818 :                   c(mems,nlh) = a(0,m,l)
     148         818 :                   lmnu(mems,nlh) = lm0
     149             :                ENDIF
     150       29084 :                DO mp=1,l
     151       15656 :                   IF( abs(a(mp,m,l)).GT.del) THEN
     152        1666 :                      mems = mems + 1
     153        1666 :                      IF(mems>memd)  CALL juDFT_error("memd too small"
     154           0 :      +                    ,calledby ="lhcal")
     155        1666 :                      c(mems,nlh) = a(mp,m,l)
     156        1666 :                      lmnu(mems,nlh) = lm0 + mp
     157        1666 :                      mems = mems + 1
     158        1666 :                      IF (mems>memd)  CALL juDFT_error("memd too small"
     159           0 :      +                    ,calledby ="lhcal")
     160        1666 :                      c(mems,nlh) = ((-1.)**mp)*conjg(c(mems-1,nlh))
     161        1666 :                      lmnu(mems,nlh) = lm0 - mp
     162             :                   ENDIF
     163             :                ENDDO
     164        2228 :                mem(nlh) = mems
     165             :             ELSE
     166       77530 :                a(0:l,m,l) = cmplx(0.0,0.0)
     167             :             ENDIF
     168             :          ENDDO   ! m = 0, l2
     169             :       ENDDO      ! l = 0, lmax
     170             : !
     171             : !===>    test of lattice harmonics using an arbitary point
     172             : !
     173         108 :       v(1) = sqrt(2.0)
     174         108 :       v(2) = sqrt(5.0)
     175         108 :       v(3) = sqrt(17.0)
     176             : !---> generate lattice harmonic for this point
     177             :       CALL ylm4(
     178             :      >          lmax,v,
     179         108 :      <          ylm)
     180        2336 :       DO lh = 1, nlh
     181        2228 :          ylms(lh) = cmplx(0.0,0.0)
     182        6486 :          DO m = 1, mem(lh)
     183        6378 :             ylms(lh) = ylms(lh) + c(m,lh)*ylm(lmnu(m,lh))
     184             :          ENDDO
     185             :       ENDDO
     186             : !---> rotate point and generate lattice harmonic
     187        1044 :       DO n = 2, nrot
     188         936 :          vr = matmul( orth(:,:,n), v )
     189             :          CALL ylm4(
     190             :      >             lmax,vr,
     191         936 :      <             ylm)
     192       14636 :          DO lh = 1, nlh
     193       13592 :             ylmr(lh) = cmplx(0.0,0.0)
     194       37110 :             DO m = 1, mem(lh)
     195       37110 :                ylmr(lh) = ylmr(lh) + c(m,lh)*ylm(lmnu(m,lh))
     196             :             ENDDO
     197       14528 :             IF ( abs(ylms(lh)-ylmr(lh)).GT.del ) THEN
     198           0 :                WRITE (6,'(/," error for operation",i3)') n
     199           0 :                WRITE (6,'(  " lattice harmonic   ",i3)') lh
     200           0 :                WRITE (6,'(4f12.6)') ylms(lh),ylmr(lh)
     201           0 :                 CALL juDFT_error("k_lv(Rr)",calledby="lhcal")
     202             :             ENDIF
     203             :          ENDDO
     204             :       ENDDO
     205             : !$OMP END MASTER
     206         108 :       RETURN
     207             :       END SUBROUTINE lhcal
     208             :       END MODULE m_lhcal

Generated by: LCOV version 1.13