checkdop.F90 9.04 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
      MODULE m_checkdop
      CONTAINS
        SUBROUTINE checkdop(&
             &                    p,np,n,na,ivac,iflag,jsp,cdn,&
             &                    DIMENSION,atoms,sphhar,stars,sym,&
             &                    vacuum,cell,oneD,&
             &                    fpw,fr,fxy,fz)
          ! ************************************************************
          !     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
          ! ************************************************************

          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 ..
          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

          INTEGER, INTENT (IN) :: iflag,ivac,n,na,np,jsp    
          LOGICAL, INTENT (IN) :: cdn 
          !-odim
          !+odim
          !     .. Array Arguments ..
          REAL,    INTENT (IN) :: p(3,DIMENSION%nspd)  
Daniel Wortmann's avatar
Daniel Wortmann committed
41
          REAL,    INTENT (IN) :: fr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd),fz(vacuum%nmzd,2,DIMENSION%jspd)
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 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 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
          COMPLEX, INTENT (IN) :: fpw(stars%n3d,DIMENSION%jspd),fxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,DIMENSION%jspd)
          !     ..
          !     .. Local Scalars ..
          REAL av,dms,rms,s,ir2,help,phi
          INTEGER i,j,k,lh,mem,nd,lm,ll1,nopa ,gz,m
          COMPLEX ic
          !     ..
          !     .. Local Arrays ..
          COMPLEX sf2(stars%n2d),sf3(stars%n3d),ylm( (atoms%lmaxd+1)**2 )
          REAL rcc(3),v1(DIMENSION%nspd),v2(DIMENSION%nspd),x(3),ri(3)
          !     ..
          !     ..
#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
                   v1(j) = v1(j) + REAL(fpw(k,jsp)*sf3(k))*stars%nstr(k)
                ENDDO
             ENDDO
             !     ---> vacuum part
             IF (cdn) THEN
                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
                   v2(j) = fz(1,ivac,jsp)
                   DO  k = 2,stars%ng2
                      v2(j) = v2(j) + REAL(fxy(1,k-1,ivac,jsp)*sf2(k))*stars%nstr2(k)
                   ENDDO
                ELSE
                   !-odim
                   v2(j) = fz(1,ivac,jsp)
                   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)
                      v2(j) = v2(j) + REAL(fxy(1,k-1,ivac,jsp)*&
                           &           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
123
             CALL fitchk(v1(:np),v2(:np),av,rms,dms)
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 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 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
             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
                v1(j) = v1(j) + REAL(fpw(k,jsp)*sf3(k))*stars%nstr(k)
             ENDDO
          ENDDO
          !     ----> m.t. part
          IF (cdn) THEN
             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
          IF (cdn) ir2 = 1.0 / ( atoms%rmt(n)*atoms%rmt(n) )
          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
                help = help + fr(atoms%jri(n),lh,n,jsp) * s
             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