MODULE m_stepf USE m_juDFT USE m_cdn_io CONTAINS SUBROUTINE stepf(sym,stars,atoms,oneD, input,cell, vacuum,mpi) ! !********************************************************************* ! 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 ! !********************************************************************* #include"cpp_double.h" 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(INOUT) :: atoms TYPE(t_oneD),INTENT(IN) :: oneD TYPE(t_input),INTENT(IN) :: input TYPE(t_cell),INTENT(INOUT) :: cell TYPE(t_vacuum),INTENT(IN) :: vacuum TYPE(t_mpi),INTENT(IN) :: mpi ! .. ! .. 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 INTEGER i,k,n,na,nn,i1,i2,i3,ic,ifft2d,ifftd,kk INTEGER ic1,ic2,ic3,icc,im1,im2,im3,loopstart LOGICAL l_error ! .. ! .. Local Arrays .. COMPLEX sf(stars%ng3) REAL g(3),gm(3),fJ REAL, ALLOCATABLE :: bfft(:) INTEGER, ALLOCATABLE :: icm(:,:,:) 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(:) 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 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 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 gm(3)=REAL(i3) 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 IF (i3 == 0 .AND. i2 == 0 ) THEN loopstart=1 ELSE loopstart=0 ENDIF DO i1=loopstart,3*stars%mx1-1 ic=i1 + 3*stars%mx1*i2 + 9*stars%mx1*stars%mx2*i3 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)) #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 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 #ifdef CPP_MPI ufft_local(ic) = help * r_c #else stars%ufft(ic) = help * r_c #endif 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 #ifdef CPP_MPI ufft_local(ic) = help * REAL(c_c) bfft_local(ic) = help * AIMAG(c_c) #else stars%ufft(ic) = help * REAL(c_c) bfft(ic) = help * AIMAG(c_c) #endif ENDIF IF (((i3.EQ.3*stars%mx3/2).OR. (i2.EQ.3*stars%mx2/2)).OR. (i1.EQ.3*stars%mx1/2)) THEN #ifdef CPP_MPI ufft_local(ic)=0.0 bfft_local(ic)=0.0 #else stars%ufft(ic)=0.0 bfft(ic)=0.0 #endif ENDIF !-odim IF (oneD%odd%d1) THEN !!!! MPI version is not tested yet !!!! IF (ic.LT.9*stars%mx1*stars%mx2 .AND. ic.NE.0) THEN 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) #ifdef CPP_MPI ufft_local(ic) = ufft_local(ic) +2*cell%vol*fJ/(gr*cell%z1*cell%omtil) #else stars%ufft(ic) = stars%ufft(ic) +2*cell%vol*fJ/(gr*cell%z1*cell%omtil) #endif END IF END IF !+odim 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) ic=ic+1 ENDDO ENDDO ENDDO ! ! --> add film-contributions ! IF (input%film .AND. .NOT.oneD%odd%d1) THEN ifft2d=9*stars%mx1*stars%mx2 stars%ufft(0)=stars%ufft(0)+cell%vol*inv_omtil-1.0 DO i3=1,3*stars%mx3-1 gm(3)=REAL(i3) IF ( gm(3) > 1.5*stars%mx3 ) gm(3)=gm(3)-3.0*stars%mx3 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 ELSEIF (oneD%odd%d1) THEN !-odim stars%ufft(0) = stars%ufft(0)+cell%vol*inv_omtil-1.0 !+odim ENDIF ! ! --> make fft ! IF (sym%invs) bfft=0.0 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) DEALLOCATE ( bfft , icm ) #ifdef CPP_MPI DEALLOCATE ( bfft_local, ufft_local , icm_local ) #endif CALL writeStepfunction(stars) ENDIF ! (mpi%irank == 0) END SUBROUTINE stepf END MODULE m_stepf