Skip to content

GitLab

Commit 458282b8 authored by Marco Govoni's avatar Marco Govoni
Browse files

Introduced function3d

parent d2434cf5
Hide whitespace changes
Inline Side-by-side
FFT_kernel/fft_interpolation.f90
View file @ 458282b8
...@@ -17,11 +17,11 @@ MODULE fourier_interpolation ...@@ -17,11 +17,11 @@ MODULE fourier_interpolation
! !
CONTAINS CONTAINS
! !
SUBROUTINE get_G2R_mapping (n1, n2, n3, n, nl, ndim) SUBROUTINE get_G2R_mapping (n1, n2, n3, ng, ngx, ndim, nl)
! !
! INPUT : n1, n2, n3 = dimension of arbitrary R grid ! INPUT : n1, n2, n3 = dimension of arbitrary R grid
! n = actual number of PW ! ng, ngx = actual number of PW
! ndim - dimensions of nl, 1: (n), 2:(n,-n) ! ndim = dimensions of nl, 1: (n), 2:(n,-n)
! OUTPUT : nl = computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] ) ! OUTPUT : nl = computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] )
! !
USE kinds, ONLY : DP USE kinds, ONLY : DP
...@@ -32,8 +32,8 @@ MODULE fourier_interpolation ...@@ -32,8 +32,8 @@ MODULE fourier_interpolation
! !
! I/O ! I/O
! !
INTEGER, INTENT(IN) :: n1, n2, n3, n, ndim INTEGER, INTENT(IN) :: n1, n2, n3, ng, ngx, ndim
INTEGER, INTENT(OUT) :: nl(n,ndim) ! mapping from [1, npw] to [1, n1*n2*n3] INTEGER, INTENT(OUT) :: nl(ngx,ndim) ! mapping from [1, npw] to [1, n1*n2*n3]
! !
! Workspace ! Workspace
! !
...@@ -49,7 +49,7 @@ MODULE fourier_interpolation ...@@ -49,7 +49,7 @@ MODULE fourier_interpolation
lmax = (n3 - 1) / 2 lmax = (n3 - 1) / 2
! !
DO i_dim = 1, ndim DO i_dim = 1, ndim
DO ig = 1, n DO ig = 1, ng
! !
! Get the current Miller index ! Get the current Miller index
! !
...@@ -98,15 +98,15 @@ MODULE fourier_interpolation ...@@ -98,15 +98,15 @@ MODULE fourier_interpolation
END SUBROUTINE END SUBROUTINE
! !
! !
SUBROUTINE single_fwfft_fromArbitraryRGrid (fr, n1, n2, n3, n, nx, ndim, nl, fg, igk) SUBROUTINE single_fwfft_fromArbitraryRGrid (fr, n1, n2, n3, ng, ngx, ndim, nl, fg, igk)
! !
! Note that the interpolation needs to be called by all processors within one band group ! Note that the interpolation needs to be called by all processors within one band group
! !
! FWFFT : R ---> G ! FWFFT : R ---> G
! !
! INPUT : n1, n2, n3 = dimension of arbitrary R grid ! INPUT : n1, n2, n3 = dimension of arbitrary R grid
! n = actual number of PW ! ng = actual number of PW
! nx = leading dimendion for fg ! ngx = leading dimendion for fg
! ndim = 1, 2 ! ndim = 1, 2
! fr = ONE COMPLEX array containing ONE function in R space (note that the array is not distributed, i.e. dimension = n1*n2*n3 ) ! fr = ONE COMPLEX array containing ONE function in R space (note that the array is not distributed, i.e. dimension = n1*n2*n3 )
! nl = pre-computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] ) ! nl = pre-computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] )
...@@ -119,11 +119,11 @@ MODULE fourier_interpolation ...@@ -119,11 +119,11 @@ MODULE fourier_interpolation
! !
! I/O ! I/O
! !
INTEGER, INTENT(IN) :: n1, n2, n3, n, nx, ndim INTEGER, INTENT(IN) :: n1, n2, n3, ng, ngx, ndim
INTEGER, INTENT(IN) :: nl(n,ndim) INTEGER, INTENT(IN) :: nl(ngx,ndim)
COMPLEX(DP), INTENT(IN) :: fr(n1*n2*n3) COMPLEX(DP), INTENT(IN) :: fr(n1*n2*n3)
COMPLEX(DP), INTENT(OUT) :: fg(nx) COMPLEX(DP), INTENT(OUT) :: fg(ngx)
INTEGER, INTENT(IN), OPTIONAL :: igk(n) INTEGER, INTENT(IN), OPTIONAL :: igk(ng)
! !
! Workspace ! Workspace
! !
...@@ -142,14 +142,14 @@ MODULE fourier_interpolation ...@@ -142,14 +142,14 @@ MODULE fourier_interpolation
CALL mp_bcast( fr, 0, intra_bgrp_comm ) CALL mp_bcast( fr, 0, intra_bgrp_comm )
! !
IF( PRESENT(igk) ) THEN IF( PRESENT(igk) ) THEN
DO ig=1, n DO ig=1, ng
! !
idx = nl(igk(ig),1) idx = nl(igk(ig),1)
IF (idx > 0) fg(ig) = fr(nl(igk(ig),1)) IF (idx > 0) fg(ig) = fr(nl(igk(ig),1))
! !
ENDDO ENDDO
ELSE ELSE
DO ig=1, n DO ig=1, ng
! !
idx = nl(ig,1) idx = nl(ig,1)
IF (idx > 0) fg(ig) = fr(nl(ig,1)) IF (idx > 0) fg(ig) = fr(nl(ig,1))
...@@ -157,7 +157,7 @@ MODULE fourier_interpolation ...@@ -157,7 +157,7 @@ MODULE fourier_interpolation
ENDDO ENDDO
ENDIF ENDIF
! !
DO ig = (n+1), nx DO ig = (ng+1), ngx
! !
fg(ig) = (0.0_DP,0.0_DP) fg(ig) = (0.0_DP,0.0_DP)
! !
...@@ -166,11 +166,11 @@ MODULE fourier_interpolation ...@@ -166,11 +166,11 @@ MODULE fourier_interpolation
ENDSUBROUTINE ENDSUBROUTINE
! !
! !
SUBROUTINE single_invfft_toArbitraryRGrid (fr, n1, n2, n3, n, nx, ndim, nl, fg, igk) SUBROUTINE single_invfft_toArbitraryRGrid (fr, n1, n2, n3, ng, ngx, ndim, nl, fg, igk)
! !
! Note that the interpolation needs to be called by all processors within one band group ! Note that the interpolation needs to be called by all processors within one band group
! !
! FWFFT : R ---> G ! INVFFT : R <--- G
! !
! INPUT : n1, n2, n3 = dimension of arbitrary R grid ! INPUT : n1, n2, n3 = dimension of arbitrary R grid
! n = actual number of PW ! n = actual number of PW
...@@ -187,11 +187,11 @@ MODULE fourier_interpolation ...@@ -187,11 +187,11 @@ MODULE fourier_interpolation
! !
! I/O ! I/O
! !
INTEGER, INTENT(IN) :: n1, n2, n3, n, nx, ndim INTEGER, INTENT(IN) :: n1, n2, n3, ng, ngx, ndim
INTEGER, INTENT(IN) :: nl(n,ndim) INTEGER, INTENT(IN) :: nl(ng,ndim)
COMPLEX(DP), INTENT(OUT) :: fr(n1*n2*n3) COMPLEX(DP), INTENT(OUT) :: fr(n1*n2*n3)
COMPLEX(DP), INTENT(IN):: fg(nx) COMPLEX(DP), INTENT(IN):: fg(ngx)
INTEGER, INTENT(IN), OPTIONAL :: igk(n) INTEGER, INTENT(IN), OPTIONAL :: igk(ng)
! !
! Workspace ! Workspace
! !
...@@ -200,7 +200,7 @@ MODULE fourier_interpolation ...@@ -200,7 +200,7 @@ MODULE fourier_interpolation
fr = 0._DP fr = 0._DP
! !
IF( PRESENT(igk) ) THEN IF( PRESENT(igk) ) THEN
DO ig=1, n DO ig=1, ng
! !
idx = nl(igk(ig),1) idx = nl(igk(ig),1)
IF (idx > 0) THEN IF (idx > 0) THEN
...@@ -209,7 +209,7 @@ MODULE fourier_interpolation ...@@ -209,7 +209,7 @@ MODULE fourier_interpolation
! !
ENDDO ENDDO
ELSE ELSE
DO ig=1, n DO ig=1, ng
! !
idx = nl(ig,1) idx = nl(ig,1)
IF (idx > 0) THEN IF (idx > 0) THEN
...@@ -218,7 +218,7 @@ MODULE fourier_interpolation ...@@ -218,7 +218,7 @@ MODULE fourier_interpolation
! !
ENDDO ENDDO
IF ( ndim == 2 ) THEN IF ( ndim == 2 ) THEN
DO ig=1, n DO ig=1, ng
! !
idx = nl(ig,2) idx = nl(ig,2)
IF (idx > 0) THEN IF (idx > 0) THEN
......
IO_kernel/Makefile
View file @ 458282b8
...@@ -19,6 +19,7 @@ IO_KERNEL_OBJS = \ ...@@ -19,6 +19,7 @@ IO_KERNEL_OBJS = \
mpiio.o \ mpiio.o \
my_mkdir.o \ my_mkdir.o \
cubefile.o \ cubefile.o \
function3d.o \
pdep_io.o \ pdep_io.o \
pdep_db.o \ pdep_db.o \
wfreq_db.o \ wfreq_db.o \
......
IO_kernel/function3d.f90 0 → 100644
View file @ 458282b8
!
! Copyright (C) 2015-2017 M. Govoni
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! This file is part of WEST.
!
! Contributors to this file:
! Marco Govoni
!
! -------------------------------------------------------------------
MODULE function3d
! -----------------------------------------------------------------
!
IMPLICIT NONE
!
CONTAINS
!
!-----------------------------------------------------------------
SUBROUTINE write_function3d ( fname, nx, ny, nz, ng, ngx, funct3d_g )
! -----------------------------------------------------------------
!
USE kinds, ONLY : DP
USE cell_base, ONLY : celldm, at
USE control_flags, ONLY : gamma_only
USE mp_bands, ONLY : me_bgrp
USE base64_module
!
IMPLICIT NONE
!
! I/O
!
CHARACTER(LEN=*),INTENT(IN) :: fname
INTEGER, INTENT(IN) :: nx, ny, nz, ng, ngx
COMPLEX(DP),INTENT(IN) :: funct3d_g(ngx)
!
! Workspace
!
INTEGER :: iu, ndim, nbytes, nlen, nmaps, i
COMPLEX(DP),ALLOCATABLE :: funct3d_r_complex(:)
REAL(DP),ALLOCATABLE :: funct3d_r_double(:)
INTEGER, ALLOCATABLE :: nl(:,:)
CHARACTER(LEN=14) :: lab(3)
CHARACTER(LEN=:),ALLOCATABLE :: charbase64
CHARACTER(LEN=:),ALLOCATABLE :: ctype
!
! 1) Fourier interpolate funct3_g --> funct3d_r
!
IF( gamma_only ) THEN
nmaps = 2
ELSE
nmaps = 1
ENDIF
ALLOCATE( nl(ng,nmaps) )
CALL get_G2R_mapping (nx, ny, nz, ng, ngx, nl, nmaps)
ALLOCATE( funct3d_r_complex(nx*ny*nz) )
CALL single_invfft_toArbitraryRGrid (funct3d_r_complex, nx, ny, nz, ng, ngx, ndim, nl, funct3d_g)
DEALLOCATE( nl )
!
IF( me_bgrp == 0 ) THEN
!
! 2) Encode
!
IF( gamma_only ) THEN
ALLOCATE( funct3d_r_double( nx*ny*nz ) )
funct3d_r_double = REAL(funct3d_r_complex(:),KIND=DP)
DEALLOCATE(funct3d_r_complex)
ndim = nx*ny*nz
nbytes = SIZEOF(funct3d_r_double(1)) * ndim
nlen = lenbase64(nbytes)
ALLOCATE(CHARACTER(LEN=nlen) :: charbase64)
IF (.NOT. islittleendian()) CALL base64_byteswap_double(nbytes,funct3d_r_double(1:ndim))
CALL base64_encode_double(funct3d_r_double(1:ndim), ndim, charbase64)
DEALLOCATE(funct3d_r_double)
ctype = "double"
ELSE
ndim = nx*ny*nz
nbytes = SIZEOF(funct3d_r_complex(1)) * ndim
nlen = lenbase64(nbytes)
ALLOCATE(CHARACTER(LEN=nlen) :: charbase64)
IF (.NOT. islittleendian()) CALL base64_byteswap_complex(nbytes,funct3d_r_complex(1:ndim))
CALL base64_encode_complex(funct3d_r_complex(1:ndim), ndim, charbase64)
DEALLOCATE(funct3d_r_complex)
ctype = "complex"
ENDIF
!
! 3) Write
!
OPEN(NEWUNIT=iu,FILE=TRIM(ADJUSTL(fname)))
!
WRITE(iu,'(a)') '<?xml version="1.0" encoding="UTF-8"?>'
WRITE(iu,'(a)') '<fpmd:function3d xmlns:fpmd="http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0"'
WRITE(iu,'(a)') 'xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"'
WRITE(iu,'(a)') 'xsi:schemaLocation="http://www.quantum-simulation.org/ns/fpmd/fpmd-1.0 function3d.xsd"'
WRITE(iu,'(a)') 'name="delta_rho">'
DO i = 1, 3
WRITE(lab(i),'(f14.6)') celldm(1) * at(i,1)
ENDDO
WRITE(iu,'(a)') '<domain a="',TRIM(ADJUSTL(lab(1))), TRIM(ADJUSTL(lab(2))), TRIM(ADJUSTL(lab(3))),'"'
DO i = 1, 3
WRITE(lab(i),'(f14.6)') celldm(1) * at(i,2)
ENDDO
WRITE(iu,'(a)') 'b="',TRIM(ADJUSTL(lab(1))), TRIM(ADJUSTL(lab(2))), TRIM(ADJUSTL(lab(3))),'"'
DO i = 1, 3
WRITE(lab(i),'(f14.6)') celldm(1) * at(i,3)
ENDDO
WRITE(iu,'(a)') 'c="',TRIM(ADJUSTL(lab(1))), TRIM(ADJUSTL(lab(2))), TRIM(ADJUSTL(lab(3))),'"/>'
WRITE(lab(1),'(i14)') nx
WRITE(lab(2),'(i14)') ny
WRITE(lab(3),'(i14)') nz
WRITE(iu,'(a)') '<grid nx="',TRIM(ADJUSTL(lab(1))),'" ny="',TRIM(ADJUSTL(lab(2))),'" nz="',TRIM(ADJUSTL(lab(3))),'"/>'
WRITE(iu,'(a)') '<grid_function type="',ctype,'" nx="',TRIM(ADJUSTL(lab(1))),'" ny="',TRIM(ADJUSTL(lab(2))), &
&'" nz="',TRIM(ADJUSTL(lab(3))),' encoding="base64"/>'
CALL write_long_string(iu,charbase64)
WRITE(iu,'(a)') '</grid_function>'
WRITE(iu,'(a)') '</function3d>'
!
CLOSE(iu)
!
DEALLOCATE(charbase64)
DEALLOCATE(ctype)
!
ELSE
DEALLOCATE( funct3d_r_complex )
ENDIF
!
END SUBROUTINE
!
!-----------------------------------------------------------------
SUBROUTINE read_function3d ( fname, nx, ny, nz, ng, ngx, funct3d_g )
! -----------------------------------------------------------------
!
USE kinds, ONLY : DP
USE control_flags, ONLY : gamma_only
USE mp_bands, ONLY : me_bgrp
USE base64_module
!
IMPLICIT NONE
!
! I/O
!
CHARACTER(LEN=*),INTENT(IN) :: fname
INTEGER,INTENT(OUT) :: nx, ny, nz
INTEGER,INTENT(IN) :: ng, ngx
COMPLEX(DP),INTENT(OUT) :: funct3d_g(ngx)
!
! Workspace
!
INTEGER :: iu, i, nlen, ndim, nbytes, nmaps
COMPLEX(DP),ALLOCATABLE :: funct3d_r_complex(:)
REAL(DP),ALLOCATABLE :: funct3d_r_double(:)
INTEGER :: nlines, j, is, ie, ios, nline_start, nline_end
CHARACTER(LEN=256) :: buffline
CHARACTER(LEN=:),ALLOCATABLE :: bufftag
CHARACTER(LEN=:),ALLOCATABLE :: buff2
CHARACTER(LEN=:),ALLOCATABLE :: charbase64
LOGICAL :: lread
LOGICAL :: lstop
INTEGER, ALLOCATABLE :: nl(:,:)
CHARACTER(LEN=:),ALLOCATABLE :: ctype
!
IF( me_bgrp == 0 ) THEN
!
OPEN(NEWUNIT=iu,FILE=TRIM(ADJUSTL(fname)))
!
DO
READ(iu,'(a)',IOSTAT=ios) buffline
IF( ios /=0 ) EXIT
IF( INDEX(buffline,"<grid_function ") /= 0 ) THEN
bufftag = TRIM(buffline)
EXIT
ENDIF
ENDDO
!
DO WHILE ( INDEX(bufftag,">") == 0 )
READ(iu,'(a)',IOSTAT=ios) buffline
IF( ios /=0 ) EXIT
bufftag = bufftag // TRIM(buffline)
ENDDO
!
IF( INDEX(bufftag,">") == 0 ) THEN
lread = .false.
ELSE
is = INDEX(bufftag,'nx="') + 4
buff2 = bufftag(is:)
ie = INDEX(buff2,'"') + is - 2
READ(bufftag(is:ie),*) nx
is = INDEX(bufftag,'ny="') + 4
buff2 = bufftag(is:)
ie = INDEX(buff2,'"') + is - 2
READ(bufftag(is:ie),*) ny
is = INDEX(bufftag,'nz="') + 4
buff2 = bufftag(is:)
ie = INDEX(buff2,'"') + is - 2
READ(bufftag(is:ie),*) nz
is = INDEX(bufftag,'type="') + 6
buff2 = bufftag(is:)
ie = INDEX(buff2,'"') + is - 2
ctype = bufftag(is:ie)
lread = .true.
ENDIF
!
lstop = .FALSE.
!
IF( lread ) THEN
ndim = nx*ny*nz
SELECT CASE (ctype)
CASE("double")
nbytes = SIZEOF(1._DP) * ndim
CASE("complex")
nbytes = SIZEOF(CMPLX(1._DP,0._DP,KIND=DP)) * ndim
CASE DEFAULT
END SELECT
nlen = lenbase64(nbytes)
ALLOCATE(CHARACTER(LEN=nlen) :: charbase64)
!
DO
READ(iu,'(a)',IOSTAT=ios) buffline
IF( ios /=0 ) EXIT
IF(INDEX(buffline,"</grid_function>") /= 0) THEN
lstop = .TRUE.
EXIT
ENDIF
IF( .NOT. ALLOCATED( charbase64 )) THEN
charbase64 = TRIM( buffline )
ELSE
charbase64 = charbase64 // TRIM( buffline )
ENDIF
ENDDO
ELSE
CALL errore("","Could not start tag",1)
ENDIF
!
IF( .NOT. lstop ) THEN
CALL errore("","Could not close tag",1)
ENDIF
!
CLOSE(iu)
!
ALLOCATE( funct3d_r_complex(1:ndim) )
!
SELECT CASE(ctype)
CASE("double")
ALLOCATE( funct3d_r_double(1:ndim) )
CALL base64_decode_double(charbase64(1:nlen), ndim, funct3d_r_double(1:ndim))
DEALLOCATE( charbase64 )
IF (.NOT. islittleendian()) CALL base64_byteswap_double(nbytes,funct3d_r_double(1:ndim))
funct3d_r_complex(:) = CMPLX( funct3d_r_double(:), 0._DP, KIND=DP)
DEALLOCATE( funct3d_r_double )
CASE("complex")
CALL base64_decode_complex(charbase64(1:nlen), ndim, funct3d_r_complex(1:ndim))
DEALLOCATE( charbase64 )
IF (.NOT. islittleendian()) CALL base64_byteswap_complex(nbytes,funct3d_r_complex(1:ndim))
CASE DEFAULT
END SELECT
!
ENDIF
!
! 1) F interpolate funct3_r --> funct3d_g
!
IF( gamma_only ) THEN
nmaps = 2
ELSE
nmaps = 1
ENDIF
ALLOCATE( nl(ng,nmaps) )
CALL get_G2R_mapping (nx, ny, nz, ng, ngx, nl, nmaps)
CALL single_fwfft_toArbitraryRGrid (funct3d_r_complex, nx, ny, nz, ng, ngx, ndim, nl, funct3d_g)
DEALLOCATE( nl )
DEALLOCATE( funct3d_r_complex )
!
END SUBROUTINE
!
!
SUBROUTINE write_long_string(iu,longstring)
!
! Write a long string on multiple lines (each line has a max of 72 charachter)
! The unit "iu" is NOT opened and closed here
!
IMPLICIT NONE
!
! I/O
!
INTEGER,INTENT(IN) :: iu
CHARACTER(LEN=*) :: longstring
!
! Workspace
!
INTEGER :: j, nlines, thislen
INTEGER, PARAMETER :: maxlen = 72
!
thislen = LEN(longstring)
nlines = thislen / maxlen
IF( MOD( thislen, maxlen ) > 0 ) nlines = nlines + 1
DO j = 1, nlines
WRITE(iu,'(a)') longstring((j-1)*72+1:MIN((j)*72,maxlen))
ENDDO
!
END SUBROUTINE
!
END MODULE
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment