LCOV - code coverage report
Current view: top level - global - checkdop.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 43 86 50.0 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_checkdop
       2             :       CONTAINS
       3          26 :         SUBROUTINE checkdop(&
       4          26 :              &                    p,np,n,na,ivac,iflag,jsp,&
       5             :              &                    DIMENSION,atoms,sphhar,stars,sym,&
       6             :              &                    vacuum,cell,oneD,potden)
       7             :           ! ************************************************************
       8             :           !     subroutines checks the continuity of coulomb           *
       9             :           !     potential or valence charge density                    *
      10             :           !     across the m.t. and vacuum boundaries                  *
      11             :           !                                     c.l.fu                 *
      12             :           !     (unifies vcheck and cdnchk      g.b.                   *
      13             :           ! YM:  this routine doesn't really work in the vacuum in 1D case yet
      14             :           ! ************************************************************
      15             : 
      16             :           USE m_juDFT
      17             :           USE m_starf, ONLY : starf2,starf3
      18             :           USE m_angle
      19             :           USE m_ylm
      20             :           USE m_types
      21             :           USE m_constants
      22             :           USE m_fitchk
      23             :           IMPLICIT NONE
      24             :           !     ..
      25             :           !     .. Scalar Arguments ..
      26             :           TYPE(t_dimension),INTENT(IN) :: dimension
      27             :           type(t_sphhar),intent(in)    :: sphhar      
      28             :           TYPE(t_stars),INTENT(IN)     :: stars
      29             :           TYPE(t_atoms),INTENT(IN)     :: atoms
      30             :           TYPE(t_sym),INTENT(IN)       :: sym
      31             :           TYPE(t_vacuum),INTENT(IN)    :: vacuum
      32             :           TYPE(t_oneD),INTENT(IN)      :: oneD
      33             :           TYPE(t_cell),INTENT(IN)      :: cell
      34             :           TYPE(t_potden),INTENT(IN)    :: potden
      35             : 
      36             :           INTEGER, INTENT (IN) :: iflag,ivac,n,na,np,jsp
      37             :           !-odim
      38             :           !+odim
      39             :           !     .. Array Arguments ..
      40             :           REAL,    INTENT (IN) :: p(3,DIMENSION%nspd)
      41             :           !     ..
      42             :           !     .. Local Scalars ..
      43             :           REAL av,dms,rms,s,ir2,help,phi
      44             :           INTEGER i,j,k,lh,mem,nd,lm,ll1,nopa ,gz,m
      45             :           COMPLEX ic
      46             :           LOGICAL l_cdn
      47             :           !     ..
      48             :           !     .. Local Arrays ..
      49         104 :           COMPLEX sf2(stars%ng2),sf3(stars%ng3),ylm( (atoms%lmaxd+1)**2 )
      50          52 :           REAL rcc(3),v1(DIMENSION%nspd),v2(DIMENSION%nspd),x(3),ri(3)
      51             : 
      52          26 :           l_cdn = .FALSE. ! By default we assume that the input is a potential.
      53          26 :           IF (potden%potdenType.LE.0) CALL juDFT_error('unknown potden type', calledby='checkdop')
      54          26 :           IF (potden%potdenType.GT.1000) l_cdn = .TRUE. ! potdenTypes > 1000 are reserved for densities
      55             :           
      56             : 
      57             :           !     ..
      58             :           !     ..
      59             : #ifdef  __TOS_BGQ__
      60             :           RETURN
      61             : #endif
      62          26 :           ic = CMPLX(0.,1.)
      63          26 :           IF (.NOT.iflag.LT.0) THEN
      64             :              !     ---> Interstitial part
      65           0 :              DO j = 1,np
      66           0 :                 IF (.NOT.oneD%odi%d1) THEN
      67             :                    !CALL cotra1(p(:,j),rcc,cell%bmat)
      68           0 :                    rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
      69             :                    CALL starf3(&
      70             :                         &               sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,p(:,j),sym%invtab,&
      71           0 :                         &               sf3)!keep
      72             :                 ENDIF
      73             :                 !
      74           0 :                 IF (oneD%odi%d1) THEN
      75             :                    !CALL cotra1(p(:,j),rcc,cell%bmat)
      76           0 :                    rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
      77             :                    CALL starf3(&
      78             :                         &               sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,rcc,sym%invtab,&
      79           0 :                         &               sf3)!keep
      80             :                 ENDIF
      81           0 :                 v1(j) = 0.0
      82           0 :                 DO k = 1,stars%ng3
      83           0 :                    v1(j) = v1(j) + REAL(potden%pw(k,jsp)*sf3(k))*stars%nstr(k)
      84             :                 ENDDO
      85             :              ENDDO
      86             :              !     ---> vacuum part
      87           0 :              IF (l_cdn) THEN
      88           0 :                 WRITE (6,FMT=9000) ivac
      89             :              ELSE
      90           0 :                 WRITE (6,FMT=8000) ivac
      91             :              ENDIF
      92           0 :              DO  j = 1,np
      93           0 :                 IF (.NOT.oneD%odi%d1) THEN
      94             :                    CALL starf2(&
      95             :                         &           sym%nop2,stars%ng2,stars%kv2,sym%mrot,sym%symor,sym%tau,p(1,j),sym%invtab,&
      96           0 :                         &           sf2)!keep
      97           0 :                    v2(j) = potden%vacz(1,ivac,jsp)
      98           0 :                    DO  k = 2,stars%ng2
      99           0 :                       v2(j) = v2(j) + REAL(potden%vacxy(1,k-1,ivac,jsp)*sf2(k))*stars%nstr2(k)
     100             :                    ENDDO
     101             :                 ELSE
     102             :                    !-odim
     103           0 :                    v2(j) = potden%vacz(1,ivac,jsp)
     104           0 :                    phi = angle(p(1,j),p(2,j))
     105           0 :                    DO  k = 2,oneD%odi%nq2
     106           0 :                       m = oneD%odi%kv(2,k)
     107           0 :                       gz = oneD%odi%kv(1,k)
     108             :                       v2(j) = v2(j) + REAL(potden%vacxy(1,k-1,ivac,jsp)*&
     109           0 :                            &           EXP(ic*m*phi)*EXP(ic*cell%bmat(3,3)*gz*p(3,j)))*oneD%odi%nst2(k)
     110             :                    ENDDO
     111             :                    !+odim
     112             :                 END IF
     113           0 :                 IF (oneD%odi%d1) THEN
     114             :                    !CALL cotra1(p(:,j),rcc,cell%bmat)
     115           0 :                    rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
     116             : 
     117           0 :                    WRITE (6,FMT=8020)  rcc,(p(i,j),i=1,3),v1(j),v2(j)
     118             :                 ELSE
     119             :                    !CALL cotra0(p(1,j),rcc,cell%amat)
     120           0 :                    rcc=MATMUL(cell%amat,p(:,j))
     121           0 :                    WRITE (6,FMT=8020) (p(i,j),i=1,3),rcc,v1(j),v2(j)
     122             :                 ENDIF
     123             :              ENDDO
     124           0 :              CALL fitchk(v1(:np),v2(:np),av,rms,dms)
     125           0 :              WRITE (6,FMT=8030) av,rms,dms
     126           0 :              RETURN
     127             :           ENDIF
     128             :           !      ----> interstitial part
     129        9126 :           DO j = 1,np
     130             :              !CALL cotra1(p(1,j),rcc,cell%bmat)
     131        9100 :              rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
     132             : 
     133             :              CALL starf3(&
     134             :                   &               sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,rcc,sym%invtab,&
     135        9100 :                   &               sf3)!keep
     136             :              !
     137        9100 :              v1(j) = 0.0
     138     6060626 :              DO k = 1,stars%ng3
     139     6060600 :                 v1(j) = v1(j) + REAL(potden%pw(k,jsp)*sf3(k))*stars%nstr(k)
     140             :              ENDDO
     141             :           ENDDO
     142             :           !     ----> m.t. part
     143          26 :           IF (l_cdn) THEN
     144           6 :              WRITE (6,FMT=9010) n
     145             :           ELSE
     146          20 :              WRITE (6,FMT=8010) n
     147             :           ENDIF
     148          26 :           ir2 = 1.0
     149          26 :           IF (l_cdn) ir2 = 1.0 / ( atoms%rmt(n)*atoms%rmt(n) )
     150          26 :           nd = atoms%ntypsy(na)
     151          26 :           nopa = atoms%ngopr(na)
     152          26 :           IF (oneD%odi%d1) THEN
     153           0 :              nopa = oneD%ods%ngopr(na)
     154           0 :              nopa = oneD%ods%invtab(nopa)
     155             :           END IF
     156        9126 :           DO j = 1,np
     157       63700 :              DO i = 1,3
     158       36400 :                 x(i) = p(i,j) - atoms%pos(i,na)
     159             :              ENDDO
     160             :              ! new
     161        9100 :              IF (nopa.NE.1) THEN
     162             :                 !CALL cotra1(x,rcc,cell%bmat)  ! switch to internal units
     163           0 :                 rcc=MATMUL(cell%bmat,x)/tpi_const
     164             : 
     165           0 :                 DO i = 1,3               ! rotate into representative
     166           0 :                    ri(i) = 0.
     167           0 :                    DO k = 1,3
     168           0 :                       IF (oneD%odi%d1) THEN
     169           0 :                          ri(i) = ri(i) + oneD%ods%mrot(i,k,nopa)*rcc(k)
     170             :                       ELSE
     171           0 :                          ri(i) = ri(i) + sym%mrot(i,k,nopa)*rcc(k)
     172             :                       END IF
     173             :                    ENDDO
     174             :                 ENDDO
     175             :                 !CALL cotra0(ri,x,cell%amat)    !switch back to cartesian units
     176           0 :                 x=MATMUL(cell%amat,ri)
     177             : 
     178             :              END IF
     179             :              ! new
     180             :              CALL ylm4(&
     181             :                   &             atoms%lmax(n),x,&
     182        9100 :                   &             ylm)
     183        9100 :              help = 0.0
     184      154700 :              DO lh = 0,sphhar%nlh(nd)
     185      145600 :                 s = 0.0
     186      145600 :                 ll1 = sphhar%llh(lh,nd) * ( sphhar%llh(lh,nd) + 1 ) + 1
     187      373100 :                 DO mem = 1,sphhar%nmem(lh,nd)
     188      227500 :                    lm = ll1 + sphhar%mlh(mem,lh,nd)
     189      373100 :                    s = s + REAL( sphhar%clnu(mem,lh,nd)* ylm(lm) )
     190             :                 ENDDO
     191      154700 :                 help = help + potden%mt(atoms%jri(n),lh,n,jsp) * s
     192             :              ENDDO
     193        9100 :              v2(j) = help * ir2 
     194        9126 :              IF (j.LE.8) THEN
     195             :                 !CALL cotra1(p(1,j),rcc,cell%bmat)
     196         208 :                 rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
     197             : 
     198         208 :                 WRITE (6,FMT=8020) rcc, (p(i,j),i=1,3),v1(j),v2(j)
     199             :              END IF
     200             :           ENDDO
     201          26 :           CALL fitchk(v1(:np),v2(:np),av,rms,dms)
     202          26 :           WRITE (6,FMT=8030) av,rms,dms
     203             : 8000      FORMAT (/,'    int.-vac. boundary (potential): ivac=',i2,/,t10,&
     204             :                &       'int-coord',t36,'cart-coord',t57,' inter. ',t69,' vacuum ')
     205             : 8010      FORMAT (/,'    int.-m.t. boundary (potential): atom type=',i2,/,&
     206             :                &       t10,'int-coord',t36,'cart-coord',t57,' inter. ',t69,&
     207             :                &       '  m. t. ')
     208             : 8020      FORMAT (1x,2 (3f8.3,2x),2f12.6)
     209             : 8030      FORMAT (/,10x,'average value = ',f10.6,/,t11,'rms,dmx=',2f7.3,&
     210             :                &       ' per cent')
     211             : 9000      FORMAT (/,'    int.-vac. boundary (density): ivac=',i2,/,t10,&
     212             :                &       'int-coord',t36,'cart-coord',t57,' inter. ',t69,' vacuum ')
     213             : 9010      FORMAT (/,'    int.-m.t. boundary (density): atom type=',i2,/,t10,&
     214             :                &       'int-coord',t36,'cart-coord',t57,' inter. ',t69,'  m. t. ')
     215          26 :           RETURN
     216             :         END SUBROUTINE checkdop
     217             :       END MODULE m_checkdop

Generated by: LCOV version 1.13