stepf.F90 10.3 KB
Newer Older
1
      MODULE m_stepf
2
      USE m_juDFT
3
      USE m_cdn_io
4
      CONTAINS
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
        SUBROUTINE stepf(sym,stars,atoms,oneD, input,cell, vacuum)
          !
          !*********************************************************************
          !     calculates the fourier components of the interstitial step
          !     function for the reciprocal vectors of the star list.
          !           m. weinert  1986
          !*********************************************************************
          !
          !     also set up FFT of U(G) on a (-2G:+2G) grid for convolutions
          !
          !*********************************************************************
          USE m_cfft
          USE m_constants
          USE m_od_cylbes
          USE m_types
          IMPLICIT NONE
          !     ..
          TYPE(t_sym),INTENT(IN)        :: sym
          TYPE(t_stars),INTENT(INOUT)   :: stars
          TYPE(t_atoms),INTENT(IN)      :: atoms
          TYPE(t_oneD),INTENT(IN)       :: oneD
          TYPE(t_input),INTENT(IN)      :: input
          TYPE(t_cell),INTENT(IN)       :: cell
          TYPE(t_vacuum),INTENT(IN)     :: vacuum
          !     ..
          !     .. Local Scalars ..
          COMPLEX c_c,c_phs
          REAL c,dd,gs,th,inv_omtil,r_phs
          REAL g_rmt,g_sqr,help,g_abs,fp_omtil,r_c,gr,gx,gy
34
          INTEGER i,k,n,na,nn,i1,i2,i3,ic,ifft2d,ifftd,kk
35
          INTEGER ic1,ic2,ic3,icc,im1,im2,im3,loopstart
36
          LOGICAL l_error
37 38
          !     ..
          !     .. Local Arrays ..
39
          COMPLEX sf(stars%ng3)
40 41 42 43 44
          REAL g(3),gm(3),fJ
          REAL,    ALLOCATABLE :: bfft(:)
          INTEGER, ALLOCATABLE :: icm(:,:,:)
          !     ..
          !     ..
45
          !--->    if step function stored on disc, then just read it in
46
          !
47
          ifftd = 27*stars%mx1*stars%mx2*stars%mx3
48

49
          CALL readStepfunction(stars, atoms, cell, vacuum, l_error)
50 51 52
          IF(.NOT.l_error) THEN
             RETURN
          END IF
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
          IF (input%film) THEN
             dd = vacuum%dvac*cell%area/cell%omtil
             IF (oneD%odd%d1) dd = cell%vol/cell%omtil
          ELSE
             dd = 1.0
          END IF
          !--->    G=0 star
          c = 0.0
          DO  n = 1,atoms%ntype
             c = c + atoms%neq(n)*atoms%volmts(n)/cell%omtil
          ENDDO
          stars%ustep(1) = CMPLX(dd-c,0.0)
          !--->    G(parallel)=0  (for film)
          IF (input%film .AND. .NOT.oneD%odd%d1) THEN
             DO  k = 2,stars%ng3
                IF (stars%ig2(k).EQ.1) THEN
                   th = cell%bmat(3,3)*stars%kv3(3,k)*cell%z1
                   stars%ustep(k) = CMPLX(cell%vol*SIN(th)/th/cell%omtil,0.0)
                ELSE
                   stars%ustep(k) = CMPLX(0.0,0.0)
                END IF
             ENDDO
             !-odim
          ELSEIF (oneD%odd%d1) THEN
             DO k = 2,stars%ng3
                gr = 0.0
                IF (stars%kv3(3,k).EQ.0) THEN
                   kk = stars%ig2(k)
                   gr = stars%sk2(kk)
                   CALL od_cylbes(1,gr*cell%z1,fJ)
                   stars%ustep(k) = CMPLX(2.*dd*fJ/(gr*cell%z1),0.)
                ELSE
                   stars%ustep(k) =CMPLX(0.,0.)
                END IF

             ENDDO
             !+odim
          ELSE
             DO  k = 2,stars%ng3
                stars%ustep(k) = CMPLX(0.0,0.0)
             END DO
          END IF
          !--->    sphere contributions
          na = 0
          DO  n = 1,atoms%ntype
             c = 3.*atoms%volmts(n)/cell%omtil
             !-->     structure factors: loop over equivalent atoms
             na = na + 1
             DO  k = 2,stars%ng3
                th = -tpi_const* DOT_PRODUCT(stars%kv3(:,k),atoms%taual(:,na))
                sf(k) = CMPLX(COS(th),SIN(th))
             END DO
             DO  nn = 2,atoms%neq(n)
                na = na + 1
                DO  k = 2,stars%ng3
                   th = -tpi_const* DOT_PRODUCT(stars%kv3(:,k),atoms%taual(:,na))
                   sf(k) = sf(k) + CMPLX(COS(th),SIN(th))
                END DO
             END DO
             !--->    update step function
             DO  k = 2,stars%ng3
                gs = stars%sk3(k)*atoms%rmt(n)
                stars%ustep(k) = stars%ustep(k) - (c* (SIN(gs)/gs-COS(gs))/ (gs*gs))* sf(k)
             ENDDO
          ENDDO
          !
          ! --> set up stepfunction on fft-grid:
          !
122 123
          ALLOCATE (  bfft(0:27*stars%mx1*stars%mx2*stars%mx3-1) )
          im1=CEILING(1.5*stars%mx1); im2=CEILING(1.5*stars%mx2); im3=CEILING(1.5*stars%mx3) 
124 125 126 127 128 129 130 131 132 133 134 135 136
          ALLOCATE ( icm(-im1:im1,-im2:im2,-im3:im3) )
          icm = 0
          ic=1
          inv_omtil=1.0/cell%omtil
          fp_omtil=  -fpi_const*inv_omtil
          !DO first vector before loop
          stars%ufft(0)=0.0
          bfft(0)=0.0
          DO n=1,atoms%ntype
             stars%ufft(0)=stars%ufft(0)+atoms%neq(n)*atoms%volmts(n)
          ENDDO
          stars%ufft(0)=1.0-stars%ufft(0)*inv_omtil
          loopstart=1
137
          DO i3=0,3*stars%mx3-1
138
             gm(3)=REAL(i3)
139 140
             IF ( gm(3) > 1.5*stars%mx3 ) gm(3)=gm(3)-3.0*stars%mx3
             DO i2=0,3*stars%mx2-1
141
                gm(2)=REAL(i2)
142 143
                IF ( gm(2) > 1.5*stars%mx2 ) gm(2)=gm(2)-3.0*stars%mx2
                DO i1=loopstart,3*stars%mx1-1
144 145
                   loopstart=0 !all further loops start at i1=0
                   gm(1)=REAL(i1)
146
                   IF ( gm(1) > 1.5*stars%mx1 ) gm(1)=gm(1)-3.0*stars%mx1
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
                   !
                   !-> use inversion <-> c.c.
                   !
                   ic1 = NINT(gm(1)) ; ic2 = NINT(gm(2)) ; ic3 = NINT(gm(3))
                   IF ( gm(3) < 0.0 ) THEN  ! retreive from table icm()
                      icc = icm(-ic1,-ic2,-ic3)
                      !IF (icc.EQ.0) THEN
                      !  write(*,*) ic1,ic2,ic3,icc
                      !   CALL juDFT_error(" error in stepf! ",calledby="stepf")
                      !ENDIF
                      stars%ufft(ic) = stars%ufft(icc)
                      IF (.NOT.sym%invs) bfft(ic) = - bfft(icc)

                      ic=ic+1
                      CYCLE 
                   ELSE                         ! store number in table icm()
                      icm(ic1,ic2,ic3) = ic
                      IF (ic1 == im1) icm(-ic1,ic2,ic3) = ic
                      IF (ic2 == im2) icm(ic1,-ic2,ic3) = ic
                      IF ((ic1 == im1).AND.(ic2 == im2)) icm(-ic1,-ic2,ic3) = ic
                   ENDIF
                   g=MATMUL(TRANSPOSE(cell%bmat),gm)
                   g_sqr = DOT_PRODUCT(g,g)
                   g_abs = SQRT(g_sqr)
                   help = fp_omtil/g_sqr
                   IF (sym%invs) THEN
                      r_c = 0.0
                      !       Better no OpenMP, huge overhead! Parallel region is located in a
                      !       nested loop and is therefore created more then a billion times. 
                      !           U.Alekseeva 15.10.2015 
177
!!$OMP  PARALLEL DO PRIVATE(r_phs,nn,th,na,g_rmt,n) DEFAULT(SHARED) REDUCTION(+:r_c)
178 179 180 181 182 183 184 185 186 187
                      DO n=1,atoms%ntype
                         r_phs = 0.0
                         na=SUM(atoms%neq(:n-1))
                         DO nn=1,atoms%neq(n)
                            th=-tpi_const*DOT_PRODUCT(gm,atoms%taual(:,na+nn))
                            r_phs = r_phs + COS(th)
                         ENDDO
                         g_rmt = g_abs * atoms%rmt(n)
                         r_c=r_c+atoms%rmt(n)*(SIN(g_rmt)/g_rmt-COS(g_rmt))*r_phs
                      ENDDO
188
!!$OMP END PARALLEL DO
189 190 191
                      stars%ufft(ic) = help * r_c
                   ELSE
                      c_c=CMPLX(0.0,0.0)
192 193
!!$OMP  PARALLEL DO PRIVATE(c_phs,nn,th,na,g_rmt,n) DEFAULT(SHARED)REDUCTION(+:c_c)

194 195 196 197 198 199 200 201 202 203
                      DO n=1,atoms%ntype
                         c_phs = CMPLX(0.0,0.0)
                         na=SUM(atoms%neq(:n-1))
                         DO nn=1,atoms%neq(n)
                            th=-tpi_const*DOT_PRODUCT(gm,atoms%taual(:,na+nn))
                            c_phs = c_phs + EXP(CMPLX(0,th))
                         ENDDO
                         g_rmt = g_abs * atoms%rmt(n)
                         c_c=c_c+atoms%rmt(n)*(SIN(g_rmt)/g_rmt-COS(g_rmt))*c_phs
                      ENDDO
204
!!$OMP END PARALLEL DO
205 206 207
                      stars%ufft(ic) = help * REAL(c_c)
                      bfft(ic) = help * AIMAG(c_c)
                   ENDIF
208 209


210
                   IF (((i3.EQ.3*stars%mx3/2).OR. (i2.EQ.3*stars%mx2/2)).OR. (i1.EQ.3*stars%mx1/2)) THEN
211 212 213 214 215
                      stars%ufft(ic)=0.0 
                      bfft(ic)=0.0 
                   ENDIF
                   !-odim
                   IF (oneD%odd%d1) THEN
216
                      IF (ic.LT.9*stars%mx1*stars%mx2 .AND. ic.NE.0) THEN
217 218 219 220 221 222 223 224 225 226 227
                         gx = (cell%bmat(1,1)*gm(1) + cell%bmat(2,1)*gm(2))
                         gy = (cell%bmat(1,2)*gm(1) + cell%bmat(2,2)*gm(2))
                         gr = SQRT(gx**2 + gy**2)
                         CALL od_cylbes(1,gr*cell%z1,fJ)
                         stars%ufft(ic) = stars%ufft(ic) +2*cell%vol*fJ/(gr*cell%z1*cell%omtil)
                      END IF
                   END IF
                   !+odim
                   ic=ic+1
                ENDDO
             ENDDO
228
          ENDDO
229 230 231 232
          !
          ! --> add film-contributions
          !
          IF (input%film .AND. .NOT.oneD%odd%d1) THEN
233

234
             ifft2d=9*stars%mx1*stars%mx2
235
             stars%ufft(0)=stars%ufft(0)+cell%vol*inv_omtil-1.0
236

237
             DO i3=1,3*stars%mx3-1
238
                gm(3)=REAL(i3)
239
                IF ( gm(3) > 1.5*stars%mx3 ) gm(3)=gm(3)-3.0*stars%mx3
240 241 242
                th=cell%bmat(3,3)*gm(3)*cell%z1
                stars%ufft(i3*ifft2d)=stars%ufft(i3*ifft2d)+cell%vol*inv_omtil*SIN(th)/th
             ENDDO
243

244 245 246 247
          ELSEIF (oneD%odd%d1) THEN
             !-odim
             stars%ufft(0) = stars%ufft(0)+cell%vol*inv_omtil-1.0
             !+odim
248

249 250 251 252 253
          ENDIF
          !
          ! --> make fft
          !
          IF (sym%invs) bfft=0.0
254 255 256
          CALL cfft(stars%ufft,bfft,ifftd,3*stars%mx1,3*stars%mx1,+1)
          CALL cfft(stars%ufft,bfft,ifftd,3*stars%mx2,9*stars%mx1*stars%mx2,+1)
          CALL cfft(stars%ufft,bfft,ifftd,3*stars%mx3,ifftd,+1)
257

258
          DEALLOCATE ( bfft , icm )
259

260
          CALL writeStepfunction(stars)
261

262
        END SUBROUTINE stepf
263
      END MODULE m_stepf