checkdop.F90 9.19 KB
Newer Older
1 2 3
      MODULE m_checkdop
      CONTAINS
        SUBROUTINE checkdop(&
4
             &                    p,np,n,na,ivac,iflag,jsp,&
5
             &                    DIMENSION,atoms,sphhar,stars,sym,&
6
             &                    vacuum,cell,oneD,potden)
7 8 9 10 11 12 13 14 15
          ! ************************************************************
          !     subroutines checks the continuity of coulomb           *
          !     potential or valence charge density                    *
          !     across the m.t. and vacuum boundaries                  *
          !                                     c.l.fu                 *
          !     (unifies vcheck and cdnchk      g.b.                   *
          ! YM:  this routine doesn't really work in the vacuum in 1D case yet
          ! ************************************************************

16
          USE m_juDFT
17 18 19 20 21 22 23 24 25
          USE m_starf, ONLY : starf2,starf3
          USE m_angle
          USE m_ylm
          USE m_types
          USE m_constants
          USE m_fitchk
          IMPLICIT NONE
          !     ..
          !     .. Scalar Arguments ..
26 27 28 29 30 31 32 33 34
          TYPE(t_dimension),INTENT(IN) :: dimension
          type(t_sphhar),intent(in)    :: sphhar      
          TYPE(t_stars),INTENT(IN)     :: stars
          TYPE(t_atoms),INTENT(IN)     :: atoms
          TYPE(t_sym),INTENT(IN)       :: sym
          TYPE(t_vacuum),INTENT(IN)    :: vacuum
          TYPE(t_oneD),INTENT(IN)      :: oneD
          TYPE(t_cell),INTENT(IN)      :: cell
          TYPE(t_potden),INTENT(IN)    :: potden
35

36
          INTEGER, INTENT (IN) :: iflag,ivac,n,na,np,jsp
37 38 39
          !-odim
          !+odim
          !     .. Array Arguments ..
40
          REAL,    INTENT (IN) :: p(3,DIMENSION%nspd)
41 42 43 44 45
          !     ..
          !     .. Local Scalars ..
          REAL av,dms,rms,s,ir2,help,phi
          INTEGER i,j,k,lh,mem,nd,lm,ll1,nopa ,gz,m
          COMPLEX ic
46
          LOGICAL l_cdn
47 48
          !     ..
          !     .. Local Arrays ..
49
          COMPLEX sf2(stars%ng2),sf3(stars%ng3),ylm( (atoms%lmaxd+1)**2 )
50
          REAL rcc(3),v1(DIMENSION%nspd),v2(DIMENSION%nspd),x(3),ri(3)
51 52 53 54 55 56

          l_cdn = .FALSE. ! By default we assume that the input is a potential.
          IF (potden%potdenType.LE.0) CALL juDFT_error('unknown potden type', calledby='checkdop')
          IF (potden%potdenType.GT.1000) l_cdn = .TRUE. ! potdenTypes > 1000 are reserved for densities
          

57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
          !     ..
          !     ..
#ifdef  __TOS_BGQ__
          RETURN
#endif
          ic = CMPLX(0.,1.)
          IF (.NOT.iflag.LT.0) THEN
             !     ---> Interstitial part
             DO j = 1,np
                IF (.NOT.oneD%odi%d1) THEN
                   !CALL cotra1(p(:,j),rcc,cell%bmat)
                   rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
                   CALL starf3(&
                        &               sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,p(:,j),sym%invtab,&
                        &               sf3)!keep
                ENDIF
                !
                IF (oneD%odi%d1) THEN
                   !CALL cotra1(p(:,j),rcc,cell%bmat)
                   rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
                   CALL starf3(&
                        &               sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,rcc,sym%invtab,&
                        &               sf3)!keep
                ENDIF
                v1(j) = 0.0
                DO k = 1,stars%ng3
83
                   v1(j) = v1(j) + REAL(potden%pw(k,jsp)*sf3(k))*stars%nstr(k)
84 85 86
                ENDDO
             ENDDO
             !     ---> vacuum part
87
             IF (l_cdn) THEN
88 89 90 91 92 93 94 95 96 97 98
                WRITE (6,FMT=9000) ivac
                WRITE (16,FMT=9000) ivac
             ELSE
                WRITE (6,FMT=8000) ivac
                WRITE (16,FMT=8000) ivac
             ENDIF
             DO  j = 1,np
                IF (.NOT.oneD%odi%d1) THEN
                   CALL starf2(&
                        &           sym%nop2,stars%ng2,stars%kv2,sym%mrot,sym%symor,sym%tau,p(1,j),sym%invtab,&
                        &           sf2)!keep
99
                   v2(j) = potden%vacz(1,ivac,jsp)
100
                   DO  k = 2,stars%ng2
101
                      v2(j) = v2(j) + REAL(potden%vacxy(1,k-1,ivac,jsp)*sf2(k))*stars%nstr2(k)
102 103 104
                   ENDDO
                ELSE
                   !-odim
105
                   v2(j) = potden%vacz(1,ivac,jsp)
106 107 108 109
                   phi = angle(p(1,j),p(2,j))
                   DO  k = 2,oneD%odi%nq2
                      m = oneD%odi%kv(2,k)
                      gz = oneD%odi%kv(1,k)
110
                      v2(j) = v2(j) + REAL(potden%vacxy(1,k-1,ivac,jsp)*&
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
                           &           EXP(ic*m*phi)*EXP(ic*cell%bmat(3,3)*gz*p(3,j)))*oneD%odi%nst2(k)
                   ENDDO
                   !+odim
                END IF
                IF (oneD%odi%d1) THEN
                   !CALL cotra1(p(:,j),rcc,cell%bmat)
                   rcc=MATMUL(cell%bmat,p(:,j))/tpi_const

                   WRITE (6,FMT=8020)  rcc,(p(i,j),i=1,3),v1(j),v2(j)
                   WRITE (16,FMT=8020) rcc,(p(i,j),i=1,3),v1(j),v2(j)
                ELSE
                   !CALL cotra0(p(1,j),rcc,cell%amat)
                   rcc=MATMUL(cell%amat,p(:,j))
                   WRITE (6,FMT=8020) (p(i,j),i=1,3),rcc,v1(j),v2(j)
                   WRITE (16,FMT=8020) (p(i,j),i=1,3),rcc,v1(j),v2(j)
                ENDIF
             ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
128
             CALL fitchk(v1(:np),v2(:np),av,rms,dms)
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
             WRITE (6,FMT=8030) av,rms,dms
             WRITE (16,FMT=8030) av,rms,dms
             RETURN
          ENDIF
          !      ----> interstitial part
          DO j = 1,np
             !CALL cotra1(p(1,j),rcc,cell%bmat)
             rcc=MATMUL(cell%bmat,p(:,j))/tpi_const

             CALL starf3(&
                  &               sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,rcc,sym%invtab,&
                  &               sf3)!keep
             !
             v1(j) = 0.0
             DO k = 1,stars%ng3
144
                v1(j) = v1(j) + REAL(potden%pw(k,jsp)*sf3(k))*stars%nstr(k)
145 146 147
             ENDDO
          ENDDO
          !     ----> m.t. part
148
          IF (l_cdn) THEN
149 150 151 152 153 154 155
             WRITE (6,FMT=9010) n
             WRITE (16,FMT=9010) n
          ELSE
             WRITE (6,FMT=8010) n
             WRITE (16,FMT=8010) n
          ENDIF
          ir2 = 1.0
156
          IF (l_cdn) ir2 = 1.0 / ( atoms%rmt(n)*atoms%rmt(n) )
157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
          nd = atoms%ntypsy(na)
          nopa = atoms%ngopr(na)
          IF (oneD%odi%d1) THEN
             nopa = oneD%ods%ngopr(na)
             nopa = oneD%ods%invtab(nopa)
          END IF
          DO j = 1,np
             DO i = 1,3
                x(i) = p(i,j) - atoms%pos(i,na)
             ENDDO
             ! new
             IF (nopa.NE.1) THEN
                !CALL cotra1(x,rcc,cell%bmat)  ! switch to internal units
                rcc=MATMUL(cell%bmat,x)/tpi_const

                DO i = 1,3               ! rotate into representative
                   ri(i) = 0.
                   DO k = 1,3
                      IF (oneD%odi%d1) THEN
                         ri(i) = ri(i) + oneD%ods%mrot(i,k,nopa)*rcc(k)
                      ELSE
                         ri(i) = ri(i) + sym%mrot(i,k,nopa)*rcc(k)
                      END IF
                   ENDDO
                ENDDO
                !CALL cotra0(ri,x,cell%amat)    !switch back to cartesian units
                x=MATMUL(cell%amat,ri)

             END IF
             ! new
             CALL ylm4(&
                  &             atoms%lmax(n),x,&
                  &             ylm)
             help = 0.0
             DO lh = 0,sphhar%nlh(nd)
                s = 0.0
                ll1 = sphhar%llh(lh,nd) * ( sphhar%llh(lh,nd) + 1 ) + 1
                DO mem = 1,sphhar%nmem(lh,nd)
                   lm = ll1 + sphhar%mlh(mem,lh,nd)
                   s = s + REAL( sphhar%clnu(mem,lh,nd)* ylm(lm) )
                ENDDO
198
                help = help + potden%mt(atoms%jri(n),lh,n,jsp) * s
199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
             ENDDO
             v2(j) = help * ir2 
             IF (j.LE.8) THEN
                !CALL cotra1(p(1,j),rcc,cell%bmat)
                rcc=MATMUL(cell%bmat,p(:,j))/tpi_const

                WRITE (6,FMT=8020) rcc, (p(i,j),i=1,3),v1(j),v2(j)
                WRITE (16,FMT=8020) rcc, (p(i,j),i=1,3),v1(j),v2(j)
             END IF
          ENDDO
          CALL fitchk(v1(:np),v2(:np),av,rms,dms)
          WRITE (6,FMT=8030) av,rms,dms
          WRITE (16,FMT=8030) av,rms,dms
8000      FORMAT (/,'    int.-vac. boundary (potential): ivac=',i2,/,t10,&
               &       'int-coord',t36,'cart-coord',t57,' inter. ',t69,' vacuum ')
8010      FORMAT (/,'    int.-m.t. boundary (potential): atom type=',i2,/,&
               &       t10,'int-coord',t36,'cart-coord',t57,' inter. ',t69,&
               &       '  m. t. ')
8020      FORMAT (1x,2 (3f8.3,2x),2f12.6)
8030      FORMAT (/,10x,'average value = ',f10.6,/,t11,'rms,dmx=',2f7.3,&
               &       ' per cent')
9000      FORMAT (/,'    int.-vac. boundary (density): ivac=',i2,/,t10,&
               &       'int-coord',t36,'cart-coord',t57,' inter. ',t69,' vacuum ')
9010      FORMAT (/,'    int.-m.t. boundary (density): atom type=',i2,/,t10,&
               &       'int-coord',t36,'cart-coord',t57,' inter. ',t69,'  m. t. ')
          RETURN
        END SUBROUTINE checkdop
      END MODULE m_checkdop