stepf.F90 14.1 KB
Newer Older
1
      MODULE m_stepf
2
      USE m_juDFT
3
      USE m_cdn_io
4
      CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
5
        SUBROUTINE stepf(sym,stars,atoms,oneD, input,cell, vacuum,mpi)
6 7 8 9 10 11 12 13 14 15
          !
          !*********************************************************************
          !     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
          !
          !*********************************************************************
Daniel Wortmann's avatar
Daniel Wortmann committed
16
#include"cpp_double.h"
17 18 19 20 21 22 23 24
          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
Daniel Wortmann's avatar
Daniel Wortmann committed
25
          TYPE(t_atoms),INTENT(INOUT)   :: atoms
26 27
          TYPE(t_oneD),INTENT(IN)       :: oneD
          TYPE(t_input),INTENT(IN)      :: input
Daniel Wortmann's avatar
Daniel Wortmann committed
28
          TYPE(t_cell),INTENT(INOUT)    :: cell
29
          TYPE(t_vacuum),INTENT(IN)     :: vacuum
Daniel Wortmann's avatar
Daniel Wortmann committed
30
          TYPE(t_mpi),INTENT(IN)        :: mpi
31 32 33 34 35
          !     ..
          !     .. 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
36
          INTEGER i,k,n,na,nn,i1,i2,i3,ic,ifft2d,ifftd,kk
37
          INTEGER ic1,ic2,ic3,icc,im1,im2,im3,loopstart
38
          LOGICAL l_error
39 40
          !     ..
          !     .. Local Arrays ..
41
          COMPLEX sf(stars%ng3)
42 43 44
          REAL g(3),gm(3),fJ
          REAL,    ALLOCATABLE :: bfft(:)
          INTEGER, ALLOCATABLE :: icm(:,:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
45 46 47 48 49 50
          INTEGER :: i3_start, i3_end, chunk_size,leftover_size
#ifdef CPP_MPI
          INCLUDE 'mpif.h'
          INTEGER ierr
          INTEGER, ALLOCATABLE :: icm_local(:,:,:)
          REAL, ALLOCATABLE :: ufft_local(:), bfft_local(:)
51

Daniel Wortmann's avatar
Daniel Wortmann committed
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
        CALL MPI_BCAST(stars%mx1,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(stars%mx2,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(stars%mx3,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
#endif

        ifftd = 27*stars%mx1*stars%mx2*stars%mx3
        !     ..
        !     ..
        !--->    if step function stored on disc, then just read it in
        !
        l_error = .FALSE.
        IF (mpi%irank == 0) CALL readStepfunction(stars, atoms, cell, vacuum, l_error)
#ifdef CPP_MPI
        CALL MPI_BCAST(l_error,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
#endif
        IF(.NOT.l_error) THEN
           RETURN
        END IF

        IF (mpi%irank == 0) THEN
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 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
          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
Daniel Wortmann's avatar
Daniel Wortmann committed
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
        ENDIF ! (mpi%irank == 0)

        !
        ! --> set up stepfunction on fft-grid:
        !
#ifdef CPP_MPI
        CALL MPI_BCAST(atoms%ntype,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(atoms%nat,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(cell%omtil,1,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(cell%bmat,9,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(sym%invs,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(oneD%odd%d1,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(input%film,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(cell%z1,1,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(cell%vol,1,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
        IF(.NOT.ALLOCATED(atoms%neq)) ALLOCATE(atoms%neq(atoms%ntype))
        IF(.NOT.ALLOCATED(atoms%volmts)) ALLOCATE(atoms%volmts(atoms%ntype))
        IF(.NOT.ALLOCATED(atoms%taual)) ALLOCATE(atoms%taual(3,atoms%nat))
        IF(.NOT.ALLOCATED(atoms%rmt)) ALLOCATE(atoms%rmt(atoms%ntype))
        IF(.NOT.ALLOCATED(stars%ufft)) ALLOCATE(stars%ufft(0:27*stars%mx1*stars%mx2*stars%mx3-1))
        CALL MPI_BCAST(atoms%neq,size(atoms%neq),MPI_INTEGER,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(atoms%volmts,size(atoms%volmts),CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(atoms%taual,size(atoms%taual),CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
        CALL MPI_BCAST(atoms%rmt,size(atoms%rmt),CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
        ALLOCATE ( bfft_local(0:ifftd-1), ufft_local(0:ifftd-1) )
        bfft_local = 0.0
        ufft_local = 0.0
#endif
       
        ALLOCATE (  bfft(0:ifftd-1) )
        im1=CEILING(1.5*stars%mx1); im2=CEILING(1.5*stars%mx2); im3=CEILING(1.5*stars%mx3) 
        ALLOCATE ( icm(-im1:im1,-im2:im2,-im3:im3) )
        icm = 0
#ifdef CPP_MPI
        ALLOCATE ( icm_local(-im1:im1,-im2:im2,-im3:im3) )
        icm_local = 0
#endif

        inv_omtil=1.0/cell%omtil
        fp_omtil=  -fpi_const*inv_omtil
        !DO first vector before loop
        stars%ufft=0.0
        bfft=0.0

#ifdef CPP_MPI
        IF ( mpi%irank == 0 ) THEN
           DO n=1,atoms%ntype
              ufft_local(0)=ufft_local(0)+atoms%neq(n)*atoms%volmts(n)
           ENDDO
           ufft_local(0)=1.0-ufft_local(0)*inv_omtil
        ENDIF
#else
        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
#endif

        i3_start = 0
        i3_end = im3
#ifdef CPP_MPI
        chunk_size = im3/mpi%isize
        leftover_size = modulo(im3,mpi%isize)
        IF (mpi%irank < leftover_size ) THEN
            i3_start =  mpi%irank      * (chunk_size + 1)
            i3_end   = (mpi%irank + 1) * (chunk_size + 1) - 1
        ELSE
            i3_start = leftover_size * (chunk_size + 1) + (mpi%irank - leftover_size) * chunk_size
            i3_end = (i3_start + chunk_size) - 1
        ENDIF
#endif
        DO i3=i3_start,i3_end
210
             gm(3)=REAL(i3)
211
             DO i2=0,3*stars%mx2-1
212
                gm(2)=REAL(i2)
Daniel Wortmann's avatar
Daniel Wortmann committed
213 214 215 216 217 218
                IF ( 2*i2 > 3*stars%mx2 ) gm(2)=gm(2)-3.0*stars%mx2
                IF (i3 == 0 .AND. i2 == 0 ) THEN
                   loopstart=1 
                ELSE
                   loopstart=0
                ENDIF
219
                DO i1=loopstart,3*stars%mx1-1
Daniel Wortmann's avatar
Daniel Wortmann committed
220
                   ic=i1 + 3*stars%mx1*i2 + 9*stars%mx1*stars%mx2*i3
221
                   gm(1)=REAL(i1)
Daniel Wortmann's avatar
Daniel Wortmann committed
222
                   IF ( 2*i1 > 3*stars%mx1 ) gm(1)=gm(1)-3.0*stars%mx1
223 224 225 226
                   !
                   !-> use inversion <-> c.c.
                   !
                   ic1 = NINT(gm(1)) ; ic2 = NINT(gm(2)) ; ic3 = NINT(gm(3))
Daniel Wortmann's avatar
Daniel Wortmann committed
227 228 229 230 231 232 233 234 235 236 237
#ifdef CPP_MPI
                   icm_local(ic1,ic2,ic3) = ic
                   IF (ic1 == im1) icm_local(-ic1,ic2,ic3) = ic
                   IF (ic2 == im2) icm_local(ic1,-ic2,ic3) = ic
                   IF ((ic1 == im1).AND.(ic2 == im2)) icm_local(-ic1,-ic2,ic3) = ic
#else
                   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
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253
                   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
                      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
Daniel Wortmann's avatar
Daniel Wortmann committed
254 255 256
#ifdef CPP_MPI
                      ufft_local(ic) = help * r_c
#else
257
                      stars%ufft(ic) = help * r_c
Daniel Wortmann's avatar
Daniel Wortmann committed
258
#endif
259 260 261 262 263 264 265 266 267 268 269 270
                   ELSE
                      c_c=CMPLX(0.0,0.0)
                      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
Daniel Wortmann's avatar
Daniel Wortmann committed
271 272 273 274
#ifdef CPP_MPI
                      ufft_local(ic) = help * REAL(c_c)
                      bfft_local(ic) = help * AIMAG(c_c)
#else
275 276
                      stars%ufft(ic) = help * REAL(c_c)
                      bfft(ic) = help * AIMAG(c_c)
Daniel Wortmann's avatar
Daniel Wortmann committed
277
#endif
278
                   ENDIF
279 280


281
                   IF (((i3.EQ.3*stars%mx3/2).OR. (i2.EQ.3*stars%mx2/2)).OR. (i1.EQ.3*stars%mx1/2)) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
282 283 284 285
#ifdef CPP_MPI
                      ufft_local(ic)=0.0 
                      bfft_local(ic)=0.0 
#else
286 287
                      stars%ufft(ic)=0.0 
                      bfft(ic)=0.0 
Daniel Wortmann's avatar
Daniel Wortmann committed
288
#endif
289 290
                   ENDIF
                   !-odim
Daniel Wortmann's avatar
Daniel Wortmann committed
291
                   IF (oneD%odd%d1) THEN  !!!! MPI version is not tested yet !!!!
292
                      IF (ic.LT.9*stars%mx1*stars%mx2 .AND. ic.NE.0) THEN
293 294 295 296
                         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)
Daniel Wortmann's avatar
Daniel Wortmann committed
297 298 299
#ifdef CPP_MPI
                         ufft_local(ic) = ufft_local(ic) +2*cell%vol*fJ/(gr*cell%z1*cell%omtil)
#else
300
                         stars%ufft(ic) = stars%ufft(ic) +2*cell%vol*fJ/(gr*cell%z1*cell%omtil)
Daniel Wortmann's avatar
Daniel Wortmann committed
301
#endif
302 303 304
                      END IF
                   END IF
                   !+odim
Daniel Wortmann's avatar
Daniel Wortmann committed
305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
                ENDDO
             ENDDO
        ENDDO

#ifdef CPP_MPI
        CALL MPI_REDUCE(ufft_local,stars%ufft,ifftd,CPP_MPI_REAL, MPI_SUM,0,mpi%mpi_comm,ierr)
        CALL MPI_REDUCE(bfft_local,bfft,ifftd,CPP_MPI_REAL, MPI_SUM,0,mpi%mpi_comm,ierr)
        CALL MPI_REDUCE(icm_local,icm,size(icm),MPI_INTEGER, MPI_SUM,0,mpi%mpi_comm,ierr)
#endif
 
        IF (mpi%irank == 0) THEN
          ic = 9*stars%mx1*stars%mx2*(im3+1)
          DO i3=im3+1,3*stars%mx3-1
             gm(3)=REAL(i3)
             gm(3)=gm(3)-3.0*stars%mx3
             DO i2=0,3*stars%mx2-1
                gm(2)=REAL(i2)
                IF ( 2*i2 > 3*stars%mx2 ) gm(2)=gm(2)-3.0*stars%mx2
                DO i1=0,3*stars%mx1-1
                   gm(1)=REAL(i1)
                   IF ( 2*i1 > 3*stars%mx1 ) gm(1)=gm(1)-3.0*stars%mx1
                   !
                   !-> use inversion <-> c.c.
                   !
                   ic1 = NINT(gm(1)) ; ic2 = NINT(gm(2)) ; ic3 = NINT(gm(3))
                   icc = icm(-ic1,-ic2,-ic3)
                   stars%ufft(ic) = stars%ufft(icc)
                   IF (.NOT.sym%invs) bfft(ic) = - bfft(icc)
333 334 335
                   ic=ic+1
                ENDDO
             ENDDO
336
          ENDDO
337 338 339 340
          !
          ! --> add film-contributions
          !
          IF (input%film .AND. .NOT.oneD%odd%d1) THEN
341

342
             ifft2d=9*stars%mx1*stars%mx2
343
             stars%ufft(0)=stars%ufft(0)+cell%vol*inv_omtil-1.0
344

345
             DO i3=1,3*stars%mx3-1
346
                gm(3)=REAL(i3)
347
                IF ( gm(3) > 1.5*stars%mx3 ) gm(3)=gm(3)-3.0*stars%mx3
348 349 350
                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
351

352 353 354 355
          ELSEIF (oneD%odd%d1) THEN
             !-odim
             stars%ufft(0) = stars%ufft(0)+cell%vol*inv_omtil-1.0
             !+odim
356

357 358 359 360 361
          ENDIF
          !
          ! --> make fft
          !
          IF (sym%invs) bfft=0.0
362 363 364
          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)
365

366
          DEALLOCATE ( bfft , icm )
Daniel Wortmann's avatar
Daniel Wortmann committed
367 368 369
#ifdef CPP_MPI
          DEALLOCATE ( bfft_local, ufft_local , icm_local )
#endif
370

371
          CALL writeStepfunction(stars)
Daniel Wortmann's avatar
Daniel Wortmann committed
372 373 374
        ENDIF ! (mpi%irank == 0)
      END SUBROUTINE stepf
    END MODULE m_stepf