Functionality of IO_kernel/function3d.f90 is tested.

In IO_kernel/function3d.f90, several bugs are fixed.
In Modules/westcom.f90, an array named nlq is added to contain nl's of different q points.
In IO_kernel/pdep_db.f90, I use write_function3d subroutine to output eigenpotentials to file and I use read_function3d immediately to read eigenpotentials from files. Difference are computed.
In Wstat/davidson_diago.f90, I compute nl for each q point.

A lot of printing sentences for test purpose are in function3d.f90 and pdep_db.f90. This needs to be removed. The format of function3d.f90 is a little bit messy,which needs to be well organized later.

.gitignore

@@ -11,3 +11,4 @@
 *.out
 doc/_build/
 Modules/west_version.f90
+*.un~

IO_kernel/function3d.f90

@@ -22,15 +22,16 @@ MODULE function3d
    SUBROUTINE write_function3d ( dfft, fname, descriptor, ng, ngx, funct3d_g, nmaps, nl )
    ! -----------------------------------------------------------------
    !
-   USE kinds,           ONLY : DP
-   USE cell_base,       ONLY : celldm, at
-   USE control_flags,   ONLY : gamma_only
-   USE mp_bands,        ONLY : me_bgrp
-   USE fourier_interpolation, ONLY : set_nl
-   USE westcom,         ONLY : fftdriver
-   USE scatter_mod,           ONLY : gather_grid
-   USE fft_types,             ONLY : fft_type_descriptor
+   USE kinds,                       ONLY : DP
+   USE cell_base,                   ONLY : celldm, at
+   USE control_flags,               ONLY : gamma_only
+   USE mp_bands,                    ONLY : me_bgrp
+   USE mp_world,                    ONLY : mpime
+   USE westcom,                     ONLY : fftdriver
+   USE scatter_mod,                 ONLY : gather_grid
+   USE fft_types,                   ONLY : fft_type_descriptor
    USE base64_module
+   USE fourier_interpolation
    !
    IMPLICIT NONE
    !
@@ -57,25 +58,29 @@ MODULE function3d
    IF ( gamma_only ) THEN 
       SELECT CASE(fftdriver) 
       CASE('Wave') 
-         CALL single_interp_invfft_gamma(dfft,ng,ngx,funct3d_g,funct3d_r_complex,fftdriver,nl)
+         CALL single_interp_invfft_gamma(dfft,ng,ngx,funct3d_g,funct3d_r_complex,'Wave',nl)
       CASE('Dense')
-         CALL single_invfft_k(dfft,ng,ngx,funct3d_g,funct3d_r_complex,fftdriver,nl)
+         CALL single_interp_invfft_k(dfft,ng,ngx,funct3d_g,funct3d_r_complex,'Dense',nl)
       END SELECT
    ELSE
-      CALL single_invfft_k(dfft,ng,ngx,funct3d_g,funct3d_r_complex,fftdriver,nl)
+      CALL single_interp_invfft_k(dfft,ng,ngx,funct3d_g,funct3d_r_complex,'Wave',nl)
    ENDIF
    !
+   IF ( mpime == 1 ) WRITE(*,*) "funct3d_r_complex when writing: ", funct3d_r_complex(1:10)
+   !
    ALLOCATE( funct3d_r_complex_gathered(dfft%nr1x*dfft%nr2x*dfft%nr3x) )
+   funct3d_r_complex_gathered = 0.0_DP
    CALL gather_grid(dfft,funct3d_r_complex,funct3d_r_complex_gathered)
+   !IF ( mpime == 0 ) WRITE(*,*) "funct3d_r_complex_gathered when writing: ", funct3d_r_complex_gathered(:)
    !
-   IF( me_bgrp == 0 ) THEN
+   IF( dfft%mype == dfft%root ) THEN
       !
       ! 2) Encode 
       !
       IF( gamma_only ) THEN 
          ALLOCATE( funct3d_r_double( dfft%nr1x*dfft%nr2x*dfft%nr3x ) ) 
-         funct3d_r_double = REAL(funct3d_r_complex(:),KIND=DP)
-         DEALLOCATE(funct3d_r_complex)
+         funct3d_r_double = REAL(funct3d_r_complex_gathered(:),KIND=DP)
+         DEALLOCATE(funct3d_r_complex_gathered)
          ndim = dfft%nr1x*dfft%nr2x*dfft%nr3x
          nbytes = SIZEOF(funct3d_r_double(1)) * ndim
          nlen = lenbase64(nbytes)
@@ -86,12 +91,12 @@ MODULE function3d
          ctype = "double"
       ELSE
          ndim = dfft%nr1x*dfft%nr2x*dfft%nr3x
-         nbytes = SIZEOF(funct3d_r_complex(1)) * ndim
+         nbytes = SIZEOF(funct3d_r_complex_gathered(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)
+         IF (.NOT. islittleendian()) CALL base64_byteswap_complex(nbytes,funct3d_r_complex_gathered(1:ndim))
+         CALL base64_encode_complex(funct3d_r_complex_gathered(1:ndim), ndim, charbase64)
+         DEALLOCATE(funct3d_r_complex_gathered)
          ctype = "complex"
       ENDIF
       !
@@ -141,13 +146,14 @@ MODULE function3d
    SUBROUTINE read_function3d ( dfft, fname, ng, ngx, funct3d_g, nmaps, nl)
    ! -----------------------------------------------------------------
    !
-   USE kinds,           ONLY : DP
-   USE control_flags,   ONLY : gamma_only
-   USE mp,              ONLY : mp_bcast
-   USE mp_bands,        ONLY : me_bgrp, intra_bgrp_comm
-   USE fft_types,             ONLY : fft_type_descriptor
-   USE scatter_mod,           ONLY : scatter_grid
-   USE westcom,         ONLY : fftdriver
+   USE kinds,                       ONLY : DP
+   USE control_flags,               ONLY : gamma_only
+   USE mp,                          ONLY : mp_bcast
+   USE mp_bands,                    ONLY : me_bgrp, intra_bgrp_comm
+   USE mp_world,                    ONLY : mpime
+   USE fft_types,                   ONLY : fft_type_descriptor
+   USE scatter_mod,                 ONLY : scatter_grid
+   USE westcom,                     ONLY : fftdriver
    USE base64_module
    USE fourier_interpolation
    !
@@ -175,7 +181,7 @@ MODULE function3d
    LOGICAL :: lread 
    CHARACTER(LEN=:),ALLOCATABLE :: ctype 
    !
-   IF( me_bgrp == 0 ) THEN
+   IF( dfft%mype == dfft%root ) THEN
       !
       OPEN(NEWUNIT=iu,FILE=TRIM(ADJUSTL(fname)))
       !
@@ -237,8 +243,20 @@ MODULE function3d
       !
       CLOSE(iu)
       !
-      ALLOCATE( funct3d_r_complex_gathered(1:ndim) )
-      !
+   ENDIF
+   !
+   CALL mp_bcast(ndim,0,intra_bgrp_comm)
+   CALL mp_bcast(nx,0,intra_bgrp_comm)
+   CALL mp_bcast(ny,0,intra_bgrp_comm)
+   CALL mp_bcast(nz,0,intra_bgrp_comm)
+   !
+   IF( nx /= dfft%nr1x ) CALL errore('read','Wrong nx',1) 
+   IF( ny /= dfft%nr2x ) CALL errore('read','Wrong ny',1) 
+   IF( nz /= dfft%nr3x ) CALL errore('read','Wrong nz',1) 
+   !
+   ALLOCATE( funct3d_r_complex_gathered(1:ndim) )
+   !
+   IF ( dfft%mype == dfft%root ) THEN
       SELECT CASE(ctype) 
       CASE("double")
          ALLOCATE( funct3d_r_double(1:ndim) )
@@ -253,33 +271,32 @@ MODULE function3d
          IF (.NOT. islittleendian()) CALL base64_byteswap_complex(nbytes,funct3d_r_complex_gathered(1:ndim)) 
       CASE DEFAULT
       END SELECT
-      !
-   ENDIF 
-   !
-   CALL mp_bcast(ndim,0,intra_bgrp_comm)  
-   CALL mp_bcast(nx,0,intra_bgrp_comm)  
-   CALL mp_bcast(ny,0,intra_bgrp_comm)  
-   CALL mp_bcast(nz,0,intra_bgrp_comm) 
+   ENDIF
    !
-   IF( nx /= dfft%nr1x ) CALL errore('read','Wrong nx',1) 
-   IF( ny /= dfft%nr2x ) CALL errore('read','Wrong ny',1) 
-   IF( nz /= dfft%nr3x ) CALL errore('read','Wrong nz',1) 
+   ALLOCATE( funct3d_r_complex(dfft%nnr) )
    !
-   IF( .NOT. ALLOCATED(funct3d_r_complex)) ALLOCATE( funct3d_r_complex(dfft%nnr) )
+   !IF ( me_bgrp == 0 ) WRITE(*,*) fname,"size of gathered : ", SIZE(funct3d_r_complex_gathered)
+   !IF ( me_bgrp == 0 ) WRITE(*,*) fname,"size of distributed : ", SIZE(funct3d_r_complex)
    !
+   !IF ( mpime == 0 ) WRITE(*,*) "funct3d_r_complex_gathered when reading: ", funct3d_r_complex_gathered(:)
    CALL scatter_grid(dfft,funct3d_r_complex_gathered,funct3d_r_complex)
    !
+   !IF ( mpime == 1 ) WRITE(*,*) "funct3d_r_complex: ",  funct3d_r_complex(:)
+   !
    IF ( gamma_only ) THEN
       SELECT CASE(fftdriver)
       CASE('Wave')
          CALL single_interp_fwfft_gamma(dfft,ng,ngx,funct3d_r_complex,funct3d_g,fftdriver,nl)
       CASE('Dense')
-         CALL single_invfft_k(dfft,ng,ngx,funct3d_r_complex,funct3d_g,fftdriver,nl)
+         CALL single_interp_fwfft_k(dfft,ng,ngx,funct3d_r_complex,funct3d_g,fftdriver,nl)
       END SELECT
    ELSE
-      CALL single_invfft_k(dfft,ng,ngx,funct3d_r_complex,funct3d_g,fftdriver,nl)
-   ENDIF 
+      CALL single_interp_fwfft_k(dfft,ng,ngx,funct3d_r_complex,funct3d_g,fftdriver,nl)
+   ENDIF
+   !
+   !IF ( mpime == 1 ) WRITE(*,*) "funct3d_g: ",  funct3d_g(1:10)
    !
+   DEALLOCATE( funct3d_r_complex_gathered )
    DEALLOCATE( funct3d_r_complex )
    !
  END SUBROUTINE

IO_kernel/pdep_db.f90

@@ -43,7 +43,7 @@ MODULE pdep_db
        !
        WRITE(label_j,'(i9.9)') j
        WRITE(label_q,'(i9.9)') iq_
-       fname = "Q"//TRIM(ADJUSTL(label_q))//"E"//TRIM(ADJUSTL(label_j))//".json" 
+       fname = "Q"//TRIM(ADJUSTL(label_q))//"E"//TRIM(ADJUSTL(label_j))//".xml" 
        !
     END SUBROUTINE 
     !
@@ -63,7 +63,7 @@ MODULE pdep_db
       USE westcom,              ONLY : wstat_calculation,n_pdep_times,n_pdep_eigen,n_pdep_maxiter,n_dfpt_maxiter, &
                                      & n_steps_write_restart,n_pdep_restart_from_itr,n_pdep_read_from_file,trev_pdep, &
                                      & tr2_dfpt,l_deflate,l_kinetic_only,ev,dvg,west_prefix,trev_pdep_rel, &
-                                     & l_minimize_exx_if_active,l_use_ecutrho,wstat_save_dir,logfile 
+                                     & l_minimize_exx_if_active,l_use_ecutrho,wstat_save_dir,logfile,nlq,dfft_io,npwqx,npwq
       USE pdep_io,              ONLY : pdep_merge_and_write_G 
       USE io_push,              ONLY : io_push_bar
       USE distribution_center,  ONLY : pert
@@ -72,6 +72,7 @@ MODULE pdep_db
       USE cell_base,            ONLY : celldm,at,bg,tpiba
       USE gvect,                ONLY : ecutrho
       USE gvecw,                ONLY : ecutwfc
+      USE function3d,           ONLY : write_function3d,read_function3d
       !
       !
       IMPLICIT NONE
@@ -108,6 +109,11 @@ MODULE pdep_db
       CHARACTER(LEN=:),ALLOCATABLE  :: eigenpot_filename(:)
       CHARACTER(LEN=:),ALLOCATABLE  :: fname
       LOGICAL :: lexists
+      COMPLEX(DP) :: tot
+      INTEGER :: i
+      COMPLEX(DP),ALLOCATABLE :: tmp_dvg(:)
+      ALLOCATE(tmp_dvg(npwqx))
+      tmp_dvg = CMPLX(0.0_DP,0.0_DP,KIND=DP)
       !
       ! Assign defaut to optional parameters 
       !
@@ -206,8 +212,16 @@ MODULE pdep_db
          global_j = pert%l2g(local_j)
          IF(global_j>n_pdep_eigen) CYCLE
          ! 
-         fname = TRIM(ADJUSTL(wstat_save_dir)) // "/"// TRIM(ADJUSTL(eigenpot_filename(global_j))) 
-         CALL pdep_merge_and_write_G(fname,dvg(:,local_j),iq_)
+         fname = TRIM(ADJUSTL(wstat_save_dir)) // "/"// TRIM(ADJUSTL(eigenpot_filename(global_j)))
+         !CALL pdep_merge_and_write_G(fname,dvg(:,local_j),iq_)
+         CALL write_function3d ( dfft_io, fname, 'Descriptor', npwq, npwqx, dvg(:,local_j), 1, nlq )
+         CALL mp_barrier(world_comm)
+         CALL read_function3d ( dfft_io, fname, npwq, npwqx, tmp_dvg(:), 1, nlq)
+         tot = CMPLX(0.0_DP,0.0_DP)
+         DO i = 1, npwq
+            tot = tot + ABS(tmp_dvg(i) - dvg(i,local_j))
+         ENDDO
+         WRITE(*,*) "difference : ", tot, "tmp_dvg(1): ", tmp_dvg(1)
          !
       ENDDO
       !
@@ -243,7 +257,7 @@ MODULE pdep_db
     SUBROUTINE pdep_db_read( nglob_to_be_read, iq, lprintinfo )
       !------------------------------------------------------------------------
       !
-      USE westcom,             ONLY : n_pdep_eigen,ev,dvg,west_prefix,npwqx,wstat_save_dir
+      USE westcom,             ONLY : n_pdep_eigen,ev,dvg,west_prefix,npwqx,wstat_save_dir,npwq,npwqx,dfft_io,nlq
       USE io_global,           ONLY : stdout 
       USE mp,                  ONLY : mp_bcast,mp_barrier
       USE mp_world,            ONLY : world_comm,mpime,root

Modules/westcom.f90

@@ -71,6 +71,7 @@ MODULE scratch_area
   !INTEGER :: io_comm ! communicator for head of images (me_bgrp==0)
   !
   TYPE ( fft_type_descriptor ) :: dfft_io
+  INTEGER,ALLOCATABLE :: nlq(:,:) ! nl of q
   !
 END MODULE
 !

Wstat/davidson_diago.f90

@@ -461,7 +461,7 @@ SUBROUTINE davidson_diago_k ( )
   USE westcom,              ONLY : dvg,dng,n_pdep_eigen,trev_pdep,n_pdep_maxiter,n_pdep_basis,wstat_calculation,ev,conv,&
                                    & n_pdep_restart_from_itr,n_pdep_read_from_file,n_steps_write_restart,n_pdep_times,&
                                    & trev_pdep_rel,tr2_dfpt,l_is_wstat_converged, &
-                                   & ngq,npwq,igq_q,npwqx
+                                   & ngq,npwq,igq_q,npwqx,dfft_io,nlq
   USE pdep_db,              ONLY : pdep_db_write,pdep_db_read
   USE wstat_restart,        ONLY : wstat_restart_write, wstat_restart_clear, wstat_restart_read
   USE mp_world,             ONLY : mpime
@@ -473,6 +473,7 @@ SUBROUTINE davidson_diago_k ( )
                                    & update_with_vr_distr,refresh_with_vr_distr 
   USE types_bz_grid,        ONLY : q_grid
   USE types_coulomb,        ONLY : pot3D
+  USE fourier_interpolation,ONLY : set_nl
   USE dfpt_module
   !
   IMPLICIT NONE
@@ -579,6 +580,12 @@ SUBROUTINE davidson_diago_k ( )
      !
      npwq = ngq(iq)
      !
+     ! set nl of this q point
+     !
+     ALLOCATE( nlq(1,npwqx) )
+     nlq = 0
+     CALL set_nl (dfft_io, npwq, npwqx, 1, nlq, igq_q(:,iq))
+     !
      ! compute Coulomb potential
      !
      CALL pot3D%init('Wave',.TRUE.,'default',iq)
@@ -908,6 +915,8 @@ SUBROUTINE davidson_diago_k ( )
   !
   CALL stop_clock( 'chidiago' )
   !
+  DEALLOCATE(nlq)
+  !
   ENDDO QPOINTS_LOOP ! iq
   !
   DEALLOCATE( conv )