LCOV - code coverage report
Current view: top level - optional - outcdn.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 64 85 75.3 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_outcdn
       2             :       use m_constants
       3             :       USE m_types
       4             : !     ********************************************************
       5             : !     calculates the charge density at given point p(i=1,3)
       6             : !     ********************************************************
       7             :       CONTAINS
       8        6000 :       SUBROUTINE outcdn(&
       9             :      &                  p,n,na,iv,iflag,jsp,sliceplot,stars,&
      10             :      &                  vacuum,sphhar,atoms,sym,cell,oneD,&
      11        6000 :      &                  qpw,rhtxy,rho,rht,&
      12             :      &                  xdnout)
      13             : !
      14             :       use m_constants
      15             :       USE m_angle
      16             :       USE m_starf, ONLY : starf2,starf3
      17             :       USE m_ylm
      18             :       IMPLICIT NONE
      19             : !
      20             :       TYPE(t_sliceplot),INTENT(IN) :: sliceplot
      21             :       TYPE(t_stars),INTENT(IN)     :: stars
      22             :       TYPE(t_vacuum),INTENT(IN)    :: vacuum
      23             :       TYPE(t_sphhar),INTENT(IN)    :: sphhar
      24             :       TYPE(t_atoms),INTENT(IN)     :: atoms
      25             :       TYPE(t_sym),INTENT(IN)       :: sym
      26             :       TYPE(t_cell),INTENT(IN)      :: cell
      27             :       TYPE(t_oneD),INTENT(IN)      :: oneD
      28             : 
      29             : 
      30             : !     .. Scalar Arguments ..
      31             :       INTEGER, INTENT (IN) :: iflag,jsp,n,na,iv
      32             :       REAL,    INTENT (OUT) :: xdnout
      33             : !-odim
      34             : !+odim
      35             : !     ..
      36             : !     .. Array Arguments ..
      37             :       COMPLEX, INTENT (IN) :: qpw(:,:) !(stars%ng3,input%jspins)
      38             :       COMPLEX, INTENT (IN) :: rhtxy(:,:,:,:) !(vacuum%nmzxyd,oneD%odi%n2d-1,2,input%jspins)
      39             :       REAL,    INTENT (IN) :: rho(:,0:,:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
      40             :       REAL,    INTENT (IN) :: rht(:,:,:) !(vacuum%nmzd,2,input%jspins)
      41             :       REAL,    INTENT (INOUT) :: p(3)
      42             : !     ..
      43             : !     .. Local Scalars ..
      44             :       REAL delta,s,sx,xd1,xd2,xx1,xx2,rrr,phi
      45             :       INTEGER i,j,jp3,jr,k,lh,mem,nd,nopa,ivac,ll1,lm ,gzi,m
      46             : !     ..
      47             : !     .. Local Arrays ..
      48       12000 :       COMPLEX sf2(stars%ng2),sf3(stars%ng3),ylm((atoms%lmaxd+1)**2)
      49             :       REAL rcc(3),x(3)
      50             : 
      51             : 
      52        6000 :       ivac=iv
      53             :      
      54        6000 :       if (iflag.ne.1) THEN
      55        4688 :       if (iflag.ne.0) THEN
      56             : !     ---> interstitial part
      57             :       !CALL cotra1(p(1),rcc,cell%bmat)
      58        3488 :       rcc=matmul(cell%bmat,p)/tpi_const
      59             :       CALL starf3(&
      60             :      &            sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,rcc,sym%invtab,&
      61        3488 :      &            sf3)
      62             : !
      63        3488 :       xdnout=dot_product(real(qpw(:,jsp)*sf3(:)),stars%nstr)
      64        3488 :       RETURN
      65             : !     ---> vacuum part
      66             :       ENDIF
      67        1200 :       xdnout = 0.
      68             : !-odim
      69        1200 :       IF (oneD%odi%d1) THEN
      70           0 :          rrr = sqrt( p(1)**2 + p(2)**2 )
      71           0 :          phi = angle(p(1),p(2))
      72           0 :          jp3 = (rrr-cell%z1)/vacuum%delz
      73           0 :          delta = (rrr-cell%z1)/vacuum%delz - jp3
      74             : !*     we count 0 as point 1
      75           0 :          jp3 = jp3 + 1
      76           0 :          IF (jp3.LT.vacuum%nmz) THEN
      77             :             xdnout = rht(jp3,ivac,jsp) + delta*&
      78           0 :      &           (rht(jp3+1,ivac,jsp)-rht(jp3,ivac,jsp))
      79           0 :             IF (jp3.LT.vacuum%nmzxy) THEN
      80           0 :                xx1 = 0.
      81           0 :                xx2 = 0.
      82           0 :                DO  k = 2,oneD%odi%nq2
      83           0 :                   m = oneD%odi%kv(2,k)
      84           0 :                   gzi = oneD%odi%kv(1,k)
      85             :                   xx1 = xx1 + real(rhtxy(jp3,k-1,ivac,jsp)*&
      86             :      &                 exp(ImagUnit*m*phi)*exp(ImagUnit*gzi*cell%bmat(3,3)*p(3)))*&
      87           0 :      &                 oneD%odi%nst2(k)
      88             :                   xx2 = xx2 + real(rhtxy(jp3+1,k-1,ivac,jsp)*&
      89             :      &                 exp(ImagUnit*m*phi)*exp(ImagUnit*gzi*cell%bmat(3,3)*p(3)))*&
      90           0 :      &                 oneD%odi%nst2(k)
      91             :             ENDDO
      92           0 :                xdnout = xdnout + xx1 + delta* (xx2-xx1)
      93             :             END IF
      94             :          ELSE
      95             :             xdnout = 0.0
      96             :          END IF
      97             : 
      98             :       ELSE
      99             : !+odim      
     100        1200 :          IF (p(3).LT.0.0) THEN
     101         600 :             ivac = vacuum%nvac
     102         600 :             IF (sym%invs) THEN
     103        1800 :                p(1:2) = -p(1:2)
     104             :             END IF
     105         600 :             p(3) = abs(p(3))
     106             :          END IF
     107             :          !CALL cotra1(p,rcc,cell%bmat)
     108        1200 :          rcc=matmul(cell%bmat,p)/tpi_const
     109             :          CALL starf2(&
     110             :      &            sym%nop2,stars%ng2,stars%kv3,sym%mrot,sym%symor,sym%tau,rcc,sym%invtab,&
     111        1200 :      &            sf2)
     112             : !
     113        1200 :          jp3 = (p(3)-cell%z1)/vacuum%delz
     114        1200 :          delta = (p(3)-cell%z1)/vacuum%delz - jp3
     115             : !*     we count 0 as point 1
     116        1200 :          jp3 = jp3 + 1
     117        1200 :          IF (jp3.LT.vacuum%nmz) THEN
     118             :              xdnout = rht(jp3,ivac,jsp) + delta*&
     119        1200 :      &               (rht(jp3+1,ivac,jsp)-rht(jp3,ivac,jsp))
     120        1200 :             IF (jp3.LT.vacuum%nmzxy) THEN
     121             :                xx1 = 0.
     122             :                xx2 = 0.
     123      202800 :               DO  k = 2,stars%ng2
     124      100800 :                xx1 = xx1 + real(rhtxy(jp3,k-1,ivac,jsp)*sf2(k))*stars%nstr2(k)
     125             :                xx2 = xx2 + real(rhtxy(jp3+1,k-1,ivac,jsp)*sf2(k))*&
     126      102000 :      &               stars%nstr2(k)
     127             :    enddo
     128        1200 :               xdnout = xdnout + xx1 + delta* (xx2-xx1)
     129             :             END IF
     130             :          ELSE
     131           0 :             xdnout = 0.0
     132             :          END IF
     133             : !----> vacuum finishes
     134             :       ENDIF
     135             : 
     136             :       RETURN
     137             :       ENDIF
     138             : !     ----> m.t. part
     139             :       
     140        1312 :       nd = atoms%ntypsy(na)
     141        1312 :       nopa = atoms%ngopr(na)
     142        1312 :       IF (oneD%odi%d1) nopa = oneD%ods%ngopr(na)
     143        1312 :       sx = 0.0
     144        5248 :       DO  i = 1,3
     145        3936 :          x(i) = p(i) - atoms%pos(i,na)
     146        5248 :          sx = sx + x(i)*x(i)
     147             :    enddo
     148        1312 :       sx = sqrt(sx)
     149        1312 :       IF (nopa.NE.1) THEN
     150             : !... switch to internal units
     151             :          !CALL cotra1(x,rcc,cell%bmat)
     152         656 :          rcc=matmul(cell%bmat,x)/tpi_const
     153             : !... rotate into representative
     154        4592 :          DO  i = 1,3
     155        1968 :             p(i) = 0.
     156        8528 :             DO  j = 1,3
     157        7872 :               IF (.NOT.oneD%odi%d1) THEN
     158        5904 :                p(i) = p(i) + sym%mrot(i,j,nopa)*rcc(j)
     159             :               ELSE
     160           0 :                p(i) = p(i) + oneD%ods%mrot(i,j,nopa)*rcc(j)
     161             :               END IF
     162             :    enddo
     163             :    enddo
     164             : !... switch back to cartesian units
     165             :          !CALL cotra0(p,x,cell%amat)
     166         656 :          x=matmul(cell%amat,p)
     167             :       END IF
     168       25516 :       DO j = atoms%jri(n),2,-1
     169       25516 :          IF (sx.GE.atoms%rmsh(j,n)) EXIT
     170             :       ENDDO
     171        1312 :       jr = j
     172             :       CALL ylm4(&
     173             :      &          atoms%lmax(n),x,&
     174        1312 :      &          ylm)
     175        1312 :       xd1 = 0.0
     176        1312 :       xd2 = 0.0
     177       24928 :       DO  lh = 0, sphhar%nlh(nd)
     178       23616 :          ll1 = sphhar%llh(lh,nd) * ( sphhar%llh(lh,nd) + 1 ) + 1
     179       23616 :          s = 0.0
     180       59040 :          DO mem = 1,sphhar%nmem(lh,nd)
     181       35424 :            lm = ll1 + sphhar%mlh(mem,lh,nd)
     182       59040 :            s = s + real( sphhar%clnu(mem,lh,nd)*ylm(lm) )
     183             :          ENDDO
     184       23616 :          IF (sliceplot%plpot) THEN
     185           0 :             xd1 = xd1 + rho(jr,lh,n,jsp)*s
     186             :          ELSE
     187       23616 :             xd1 = xd1 + rho(jr,lh,n,jsp)*s/ (atoms%rmsh(jr,n)*atoms%rmsh(jr,n))
     188             :          END IF
     189       23616 :          IF (jr.EQ.atoms%jri(n)) CYCLE
     190       24928 :          IF (sliceplot%plpot) THEN
     191           0 :             xd2 = xd2 + rho(jr+1,lh,n,jsp)*s
     192             :          ELSE
     193             :             xd2 = xd2 + rho(jr+1,lh,n,jsp)*s/&
     194       23616 :      &            (atoms%rmsh(jr+1,n)*atoms%rmsh(jr+1,n))
     195             :          END IF
     196             :    ENDDO
     197        1312 :       IF (jr.EQ.atoms%jri(n)) THEN
     198           0 :          xdnout = xd1
     199             :       ELSE
     200             :          xdnout = xd1 + (xd2-xd1) *&
     201        1312 :      &                  (sx-atoms%rmsh(jr,n)) / (atoms%rmsh(jr+1,n)-atoms%rmsh(jr,n))
     202             :       END IF
     203             :  8000 FORMAT (2f10.6)
     204             : !
     205             :       RETURN
     206             :       END SUBROUTINE outcdn
     207             :       END MODULE m_outcdn

Generated by: LCOV version 1.13