LCOV - code coverage report
Current view: top level - init - stepf.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 174 189 92.1 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_stepf
       2             :       USE m_juDFT
       3             :       USE m_cdn_io
       4             :       CONTAINS
       5          76 :         SUBROUTINE stepf(sym,stars,atoms,oneD, input,cell, vacuum,mpi)
       6             :           !
       7             :           !*********************************************************************
       8             :           !     calculates the fourier components of the interstitial step
       9             :           !     function for the reciprocal vectors of the star list.
      10             :           !           m. weinert  1986
      11             :           !*********************************************************************
      12             :           !
      13             :           !     also set up FFT of U(G) on a (-2G:+2G) grid for convolutions
      14             :           !
      15             :           !*********************************************************************
      16             : #include"cpp_double.h"
      17             :           USE m_cfft
      18             :           USE m_constants
      19             :           USE m_od_cylbes
      20             :           USE m_types
      21             :           IMPLICIT NONE
      22             :           !     ..
      23             :           TYPE(t_sym),INTENT(IN)        :: sym
      24             :           TYPE(t_stars),INTENT(INOUT)   :: stars
      25             :           TYPE(t_atoms),INTENT(INOUT)   :: atoms
      26             :           TYPE(t_oneD),INTENT(IN)       :: oneD
      27             :           TYPE(t_input),INTENT(IN)      :: input
      28             :           TYPE(t_cell),INTENT(INOUT)    :: cell
      29             :           TYPE(t_vacuum),INTENT(IN)     :: vacuum
      30             :           TYPE(t_mpi),INTENT(IN)        :: mpi
      31             :           !     ..
      32             :           !     .. Local Scalars ..
      33             :           COMPLEX c_c,c_phs
      34             :           REAL c,dd,gs,th,inv_omtil,r_phs
      35             :           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          76 :           COMPLEX,ALLOCATABLE:: sf(:)
      42             :           REAL g(3),gm(3),fJ
      43          76 :           REAL,    ALLOCATABLE :: bfft(:)
      44          76 :           INTEGER, ALLOCATABLE :: icm(:,:,:)
      45             :           INTEGER :: i3_start, i3_end, chunk_size,leftover_size
      46             : #ifdef CPP_MPI
      47             :           INCLUDE 'mpif.h'
      48             :           INTEGER ierr
      49          76 :           INTEGER, ALLOCATABLE :: icm_local(:,:,:)
      50         152 :           REAL, ALLOCATABLE :: ufft_local(:), bfft_local(:)
      51             : 
      52          76 :         CALL MPI_BCAST(stars%mx1,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
      53          76 :         CALL MPI_BCAST(stars%mx2,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
      54          76 :         CALL MPI_BCAST(stars%mx3,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
      55             : #endif
      56             : 
      57          76 :         ifftd = 27*stars%mx1*stars%mx2*stars%mx3
      58             :         !     ..
      59             :         !     ..
      60             :         !--->    if step function stored on disc, then just read it in
      61             :         !
      62          76 :         l_error = .FALSE.
      63          76 :         IF (mpi%irank == 0) CALL readStepfunction(stars, atoms, cell, vacuum, l_error)
      64             : #ifdef CPP_MPI
      65          76 :         CALL MPI_BCAST(l_error,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
      66             : #endif
      67          76 :         IF(.NOT.l_error) THEN
      68          22 :            RETURN
      69             :         END IF
      70             : 
      71          54 :         IF (mpi%irank == 0) THEN
      72          27 :            ALLOCATE(sf(stars%ng3))
      73          27 :           IF (input%film) THEN
      74           4 :              dd = vacuum%dvac*cell%area/cell%omtil
      75           4 :              IF (oneD%odd%d1) dd = cell%vol/cell%omtil
      76             :           ELSE
      77             :              dd = 1.0
      78             :           END IF
      79             :           !--->    G=0 star
      80          27 :           c = 0.0
      81          87 :           DO  n = 1,atoms%ntype
      82          87 :              c = c + atoms%neq(n)*atoms%volmts(n)/cell%omtil
      83             :           ENDDO
      84          27 :           stars%ustep(1) = CMPLX(dd-c,0.0)
      85             :           !--->    G(parallel)=0  (for film)
      86          27 :           IF (input%film .AND. .NOT.oneD%odd%d1) THEN
      87       28080 :              DO  k = 2,stars%ng3
      88       14042 :                 IF (stars%ig2(k).EQ.1) THEN
      89          88 :                    th = cell%bmat(3,3)*stars%kv3(3,k)*cell%z1
      90          88 :                    stars%ustep(k) = CMPLX(cell%vol*SIN(th)/th/cell%omtil,0.0)
      91             :                 ELSE
      92       13950 :                    stars%ustep(k) = CMPLX(0.0,0.0)
      93             :                 END IF
      94             :              ENDDO
      95             :              !-odim
      96          23 :           ELSEIF (oneD%odd%d1) THEN
      97           0 :              DO k = 2,stars%ng3
      98           0 :                 gr = 0.0
      99           0 :                 IF (stars%kv3(3,k).EQ.0) THEN
     100           0 :                    kk = stars%ig2(k)
     101           0 :                    gr = stars%sk2(kk)
     102           0 :                    CALL od_cylbes(1,gr*cell%z1,fJ)
     103           0 :                    stars%ustep(k) = CMPLX(2.*dd*fJ/(gr*cell%z1),0.)
     104             :                 ELSE
     105           0 :                    stars%ustep(k) =CMPLX(0.,0.)
     106             :                 END IF
     107             : 
     108             :              ENDDO
     109             :              !+odim
     110             :           ELSE
     111      157687 :              DO  k = 2,stars%ng3
     112          23 :                 stars%ustep(k) = CMPLX(0.0,0.0)
     113             :              END DO
     114             :           END IF
     115             :           !--->    sphere contributions
     116          27 :           na = 0
     117          87 :           DO  n = 1,atoms%ntype
     118          60 :              c = 3.*atoms%volmts(n)/cell%omtil
     119             :              !-->     structure factors: loop over equivalent atoms
     120          60 :              na = na + 1
     121      133692 :              DO  k = 2,stars%ng3
     122      133632 :                 th = -tpi_const* DOT_PRODUCT(stars%kv3(:,k),atoms%taual(:,na))
     123      133692 :                 sf(k) = CMPLX(COS(th),SIN(th))
     124             :              END DO
     125          72 :              DO  nn = 2,atoms%neq(n)
     126          12 :                 na = na + 1
     127       19873 :                 DO  k = 2,stars%ng3
     128       19801 :                    th = -tpi_const* DOT_PRODUCT(stars%kv3(:,k),atoms%taual(:,na))
     129       19813 :                    sf(k) = sf(k) + CMPLX(COS(th),SIN(th))
     130             :                 END DO
     131             :              END DO
     132             :              !--->    update step function
     133      267351 :              DO  k = 2,stars%ng3
     134      133632 :                 gs = stars%sk3(k)*atoms%rmt(n)
     135      133692 :                 stars%ustep(k) = stars%ustep(k) - (c* (SIN(gs)/gs-COS(gs))/ (gs*gs))* sf(k)
     136             :              ENDDO
     137             :           ENDDO
     138             :         ENDIF ! (mpi%irank == 0)
     139             : 
     140             :         !
     141             :         ! --> set up stepfunction on fft-grid:
     142             :         !
     143             : #ifdef CPP_MPI
     144          54 :         CALL MPI_BCAST(atoms%ntype,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
     145          54 :         CALL MPI_BCAST(atoms%nat,1,MPI_INTEGER,0,mpi%mpi_comm,ierr)
     146          54 :         CALL MPI_BCAST(cell%omtil,1,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
     147          54 :         CALL MPI_BCAST(cell%bmat,9,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
     148          54 :         CALL MPI_BCAST(sym%invs,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
     149          54 :         CALL MPI_BCAST(oneD%odd%d1,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
     150          54 :         CALL MPI_BCAST(input%film,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
     151          54 :         CALL MPI_BCAST(cell%z1,1,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
     152          54 :         CALL MPI_BCAST(cell%vol,1,CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
     153          54 :         IF(.NOT.ALLOCATED(atoms%neq)) ALLOCATE(atoms%neq(atoms%ntype))
     154          54 :         IF(.NOT.ALLOCATED(atoms%volmts)) ALLOCATE(atoms%volmts(atoms%ntype))
     155          54 :         IF(.NOT.ALLOCATED(atoms%taual)) ALLOCATE(atoms%taual(3,atoms%nat))
     156          54 :         IF(.NOT.ALLOCATED(atoms%rmt)) ALLOCATE(atoms%rmt(atoms%ntype))
     157          54 :         IF(.NOT.ALLOCATED(stars%ufft)) ALLOCATE(stars%ufft(0:27*stars%mx1*stars%mx2*stars%mx3-1))
     158          54 :         CALL MPI_BCAST(atoms%neq,size(atoms%neq),MPI_INTEGER,0,mpi%mpi_comm,ierr)
     159          54 :         CALL MPI_BCAST(atoms%volmts,size(atoms%volmts),CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
     160          54 :         CALL MPI_BCAST(atoms%taual,size(atoms%taual),CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
     161          54 :         CALL MPI_BCAST(atoms%rmt,size(atoms%rmt),CPP_MPI_REAL,0,mpi%mpi_comm,ierr)
     162          54 :         ALLOCATE ( bfft_local(0:ifftd-1), ufft_local(0:ifftd-1) )
     163     3321486 :         bfft_local = 0.0
     164     3321486 :         ufft_local = 0.0
     165             : #endif
     166             :        
     167          54 :         ALLOCATE (  bfft(0:ifftd-1) )
     168          54 :         im1=CEILING(1.5*stars%mx1); im2=CEILING(1.5*stars%mx2); im3=CEILING(1.5*stars%mx3) 
     169          54 :         ALLOCATE ( icm(-im1:im1,-im2:im2,-im3:im3) )
     170        2792 :         icm = 0
     171             : #ifdef CPP_MPI
     172          54 :         ALLOCATE ( icm_local(-im1:im1,-im2:im2,-im3:im3) )
     173        2792 :         icm_local = 0
     174             : #endif
     175             : 
     176          54 :         inv_omtil=1.0/cell%omtil
     177          54 :         fp_omtil=  -fpi_const*inv_omtil
     178             :         !DO first vector before loop
     179          54 :         stars%ufft=0.0
     180     3321486 :         bfft=0.0
     181             : 
     182             : #ifdef CPP_MPI
     183          54 :         IF ( mpi%irank == 0 ) THEN
     184          87 :            DO n=1,atoms%ntype
     185          87 :               ufft_local(0)=ufft_local(0)+atoms%neq(n)*atoms%volmts(n)
     186             :            ENDDO
     187          27 :            ufft_local(0)=1.0-ufft_local(0)*inv_omtil
     188             :         ENDIF
     189             : #else
     190             :         DO n=1,atoms%ntype
     191             :            stars%ufft(0)=stars%ufft(0)+atoms%neq(n)*atoms%volmts(n)
     192             :         ENDDO
     193             :         stars%ufft(0)=1.0-stars%ufft(0)*inv_omtil
     194             : #endif
     195             : 
     196          54 :         i3_start = 0
     197          54 :         i3_end = im3
     198             : #ifdef CPP_MPI
     199          54 :         chunk_size = im3/mpi%isize
     200          54 :         leftover_size = modulo(im3,mpi%isize)
     201          54 :         IF (mpi%irank < leftover_size ) THEN
     202          15 :             i3_start =  mpi%irank      * (chunk_size + 1)
     203          15 :             i3_end   = (mpi%irank + 1) * (chunk_size + 1) - 1
     204             :         ELSE
     205          39 :             i3_start = leftover_size * (chunk_size + 1) + (mpi%irank - leftover_size) * chunk_size
     206          39 :             i3_end = (i3_start + chunk_size) - 1
     207             :         ENDIF
     208             : #endif
     209         725 :         DO i3=i3_start,i3_end
     210         671 :              gm(3)=REAL(i3)
     211       23153 :              DO i2=0,3*stars%mx2-1
     212       22428 :                 gm(2)=REAL(i2)
     213       22428 :                 IF ( 2*i2 > 3*stars%mx2 ) gm(2)=gm(2)-3.0*stars%mx2
     214       22428 :                 IF (i3 == 0 .AND. i2 == 0 ) THEN
     215             :                    loopstart=1 
     216             :                 ELSE
     217       22401 :                    loopstart=0
     218             :                 ENDIF
     219      861458 :                 DO i1=loopstart,3*stars%mx1-1
     220      838359 :                    ic=i1 + 3*stars%mx1*i2 + 9*stars%mx1*stars%mx2*i3
     221      838359 :                    gm(1)=REAL(i1)
     222      838359 :                    IF ( 2*i1 > 3*stars%mx1 ) gm(1)=gm(1)-3.0*stars%mx1
     223             :                    !
     224             :                    !-> use inversion <-> c.c.
     225             :                    !
     226      838359 :                    ic1 = NINT(gm(1)) ; ic2 = NINT(gm(2)) ; ic3 = NINT(gm(3))
     227             : #ifdef CPP_MPI
     228      838359 :                    icm_local(ic1,ic2,ic3) = ic
     229      838359 :                    IF (ic1 == im1) icm_local(-ic1,ic2,ic3) = ic
     230      838359 :                    IF (ic2 == im2) icm_local(ic1,-ic2,ic3) = ic
     231      838359 :                    IF ((ic1 == im1).AND.(ic2 == im2)) icm_local(-ic1,-ic2,ic3) = ic
     232             : #else
     233             :                    icm(ic1,ic2,ic3) = ic
     234             :                    IF (ic1 == im1) icm(-ic1,ic2,ic3) = ic
     235             :                    IF (ic2 == im2) icm(ic1,-ic2,ic3) = ic
     236             :                    IF ((ic1 == im1).AND.(ic2 == im2)) icm(-ic1,-ic2,ic3) = ic
     237             : #endif
     238      838359 :                    g=MATMUL(TRANSPOSE(cell%bmat),gm)
     239     3353436 :                    g_sqr = DOT_PRODUCT(g,g)
     240      838359 :                    g_abs = SQRT(g_sqr)
     241      838359 :                    help = fp_omtil/g_sqr
     242      838359 :                    IF (sym%invs) THEN
     243      442970 :                       r_c = 0.0
     244     1459887 :                       DO n=1,atoms%ntype
     245     1016917 :                          r_phs = 0.0
     246     1016917 :                          na=SUM(atoms%neq(:n-1))
     247     2477058 :                          DO nn=1,atoms%neq(n)
     248     1460141 :                             th=-tpi_const*DOT_PRODUCT(gm,atoms%taual(:,na+nn))
     249     2477058 :                             r_phs = r_phs + COS(th)
     250             :                          ENDDO
     251     1016917 :                          g_rmt = g_abs * atoms%rmt(n)
     252     1459887 :                          r_c=r_c+atoms%rmt(n)*(SIN(g_rmt)/g_rmt-COS(g_rmt))*r_phs
     253             :                       ENDDO
     254             : #ifdef CPP_MPI
     255      442970 :                       ufft_local(ic) = help * r_c
     256             : #else
     257             :                       stars%ufft(ic) = help * r_c
     258             : #endif
     259             :                    ELSE
     260      395389 :                       c_c=CMPLX(0.0,0.0)
     261     1308894 :                       DO n=1,atoms%ntype
     262      913505 :                          c_phs = CMPLX(0.0,0.0)
     263      913505 :                          na=SUM(atoms%neq(:n-1))
     264     1989294 :                          DO nn=1,atoms%neq(n)
     265     1075789 :                             th=-tpi_const*DOT_PRODUCT(gm,atoms%taual(:,na+nn))
     266     1989294 :                             c_phs = c_phs + EXP(CMPLX(0,th))
     267             :                          ENDDO
     268      913505 :                          g_rmt = g_abs * atoms%rmt(n)
     269     1308894 :                          c_c=c_c+atoms%rmt(n)*(SIN(g_rmt)/g_rmt-COS(g_rmt))*c_phs
     270             :                       ENDDO
     271             : #ifdef CPP_MPI
     272      395389 :                       ufft_local(ic) = help * REAL(c_c)
     273      395389 :                       bfft_local(ic) = help * AIMAG(c_c)
     274             : #else
     275             :                       stars%ufft(ic) = help * REAL(c_c)
     276             :                       bfft(ic) = help * AIMAG(c_c)
     277             : #endif
     278             :                    ENDIF
     279             : 
     280             : 
     281      838359 :                    IF (((i3.EQ.3*stars%mx3/2).OR. (i2.EQ.3*stars%mx2/2)).OR. (i1.EQ.3*stars%mx1/2)) THEN
     282             : #ifdef CPP_MPI
     283       59426 :                       ufft_local(ic)=0.0 
     284       59426 :                       bfft_local(ic)=0.0 
     285             : #else
     286             :                       stars%ufft(ic)=0.0 
     287             :                       bfft(ic)=0.0 
     288             : #endif
     289             :                    ENDIF
     290             :                    !-odim
     291      860787 :                    IF (oneD%odd%d1) THEN  !!!! MPI version is not tested yet !!!!
     292           0 :                       IF (ic.LT.9*stars%mx1*stars%mx2 .AND. ic.NE.0) THEN
     293           0 :                          gx = (cell%bmat(1,1)*gm(1) + cell%bmat(2,1)*gm(2))
     294           0 :                          gy = (cell%bmat(1,2)*gm(1) + cell%bmat(2,2)*gm(2))
     295           0 :                          gr = SQRT(gx**2 + gy**2)
     296           0 :                          CALL od_cylbes(1,gr*cell%z1,fJ)
     297             : #ifdef CPP_MPI
     298           0 :                          ufft_local(ic) = ufft_local(ic) +2*cell%vol*fJ/(gr*cell%z1*cell%omtil)
     299             : #else
     300             :                          stars%ufft(ic) = stars%ufft(ic) +2*cell%vol*fJ/(gr*cell%z1*cell%omtil)
     301             : #endif
     302             :                       END IF
     303             :                    END IF
     304             :                    !+odim
     305             :                 ENDDO
     306             :              ENDDO
     307             :         ENDDO
     308             : 
     309             : #ifdef CPP_MPI
     310          54 :         CALL MPI_REDUCE(ufft_local,stars%ufft,ifftd,CPP_MPI_REAL, MPI_SUM,0,mpi%mpi_comm,ierr)
     311          54 :         CALL MPI_REDUCE(bfft_local,bfft,ifftd,CPP_MPI_REAL, MPI_SUM,0,mpi%mpi_comm,ierr)
     312          54 :         CALL MPI_REDUCE(icm_local,icm,size(icm),MPI_INTEGER, MPI_SUM,0,mpi%mpi_comm,ierr)
     313             : #endif
     314             :  
     315          54 :         IF (mpi%irank == 0) THEN
     316          27 :           ic = 9*stars%mx1*stars%mx2*(im3+1)
     317         658 :           DO i3=im3+1,3*stars%mx3-1
     318         631 :              gm(3)=REAL(i3)
     319         631 :              gm(3)=gm(3)-3.0*stars%mx3
     320       21811 :              DO i2=0,3*stars%mx2-1
     321       21153 :                 gm(2)=REAL(i2)
     322       21153 :                 IF ( 2*i2 > 3*stars%mx2 ) gm(2)=gm(2)-3.0*stars%mx2
     323      813577 :                 DO i1=0,3*stars%mx1-1
     324      791793 :                    gm(1)=REAL(i1)
     325      791793 :                    IF ( 2*i1 > 3*stars%mx1 ) gm(1)=gm(1)-3.0*stars%mx1
     326             :                    !
     327             :                    !-> use inversion <-> c.c.
     328             :                    !
     329      791793 :                    ic1 = NINT(gm(1)) ; ic2 = NINT(gm(2)) ; ic3 = NINT(gm(3))
     330      791793 :                    icc = icm(-ic1,-ic2,-ic3)
     331      791793 :                    stars%ufft(ic) = stars%ufft(icc)
     332      791793 :                    IF (.NOT.sym%invs) bfft(ic) = - bfft(icc)
     333      812946 :                    ic=ic+1
     334             :                 ENDDO
     335             :              ENDDO
     336             :           ENDDO
     337             :           !
     338             :           ! --> add film-contributions
     339             :           !
     340          27 :           IF (input%film .AND. .NOT.oneD%odd%d1) THEN
     341             : 
     342           4 :              ifft2d=9*stars%mx1*stars%mx2
     343           4 :              stars%ufft(0)=stars%ufft(0)+cell%vol*inv_omtil-1.0
     344             : 
     345         276 :              DO i3=1,3*stars%mx3-1
     346         272 :                 gm(3)=REAL(i3)
     347         272 :                 IF ( gm(3) > 1.5*stars%mx3 ) gm(3)=gm(3)-3.0*stars%mx3
     348         272 :                 th=cell%bmat(3,3)*gm(3)*cell%z1
     349         276 :                 stars%ufft(i3*ifft2d)=stars%ufft(i3*ifft2d)+cell%vol*inv_omtil*SIN(th)/th
     350             :              ENDDO
     351             : 
     352          23 :           ELSEIF (oneD%odd%d1) THEN
     353             :              !-odim
     354           0 :              stars%ufft(0) = stars%ufft(0)+cell%vol*inv_omtil-1.0
     355             :              !+odim
     356             : 
     357             :           ENDIF
     358             :           !
     359             :           ! --> make fft
     360             :           !
     361          27 :           IF (sym%invs) bfft=0.0
     362          27 :           CALL cfft(stars%ufft,bfft,ifftd,3*stars%mx1,3*stars%mx1,+1)
     363          27 :           CALL cfft(stars%ufft,bfft,ifftd,3*stars%mx2,9*stars%mx1*stars%mx2,+1)
     364          27 :           CALL cfft(stars%ufft,bfft,ifftd,3*stars%mx3,ifftd,+1)
     365             : 
     366          27 :           DEALLOCATE ( bfft , icm )
     367             : #ifdef CPP_MPI
     368          27 :           DEALLOCATE ( bfft_local, ufft_local , icm_local )
     369             : #endif
     370             : 
     371          27 :           CALL writeStepfunction(stars)
     372             :         ENDIF ! (mpi%irank == 0)
     373          54 :       END SUBROUTINE stepf
     374             :     END MODULE m_stepf

Generated by: LCOV version 1.13