コンテンツに飛ぶ | ナビゲーションに飛ぶ

パーソナルツール

Navigation

現在位置: ホーム / Downloads / PHASE System Download / phase0_2019.01.01.patch

phase0_2019.01.01.patch

differences between files icon phase0_2019.01.01.patch — differences between files, 65 KB (67517 bytes)

ファイルコンテンツ

diff -rcN phase0_2019.01/src_phase/m_ES_IO.F90 phase0_2019.02/src_phase/m_ES_IO.F90
*** phase0_2019.01/src_phase/m_ES_IO.F90	2019-07-26 09:06:42.416093120 +0900
--- phase0_2019.02/src_phase/m_ES_IO.F90	2019-08-30 13:52:20.478646248 +0900
***************
*** 1905,1911 ****
  
    subroutine m_ESIO_wd_Psicoef(ipri,nfout,nf)
      integer, intent(in) :: ipri,nfout, nf
! 
      integer, parameter :: Ncol = 5
      integer :: ik, ie, ri,  nel, ig, ib,ib1,ib2,ibt,ibsize
      integer, allocatable, dimension(:) :: n_mpi
--- 1905,1911 ----
  
    subroutine m_ESIO_wd_Psicoef(ipri,nfout,nf)
      integer, intent(in) :: ipri,nfout, nf
!     integer :: ii,jj
      integer, parameter :: Ncol = 5
      integer :: ik, ie, ri,  nel, ig, ib,ib1,ib2,ibt,ibsize
      integer, allocatable, dimension(:) :: n_mpi
***************
*** 1923,1929 ****
      KPOINT: do ik = 1, kv3, af+1
         nel = min(Nw_Psicoef,iba(ik))
         allocate(wf_l(nel,kimg)); wf_l = 0.d0
!        call wd_k_points()
         e_mpi = 0.d0
         n_mpi = 0
         if(map_k(ik) == myrank_k) then
--- 1923,1935 ----
      KPOINT: do ik = 1, kv3, af+1
         nel = min(Nw_Psicoef,iba(ik))
         allocate(wf_l(nel,kimg)); wf_l = 0.d0
! !-----------------------------Modified by T.A.Ariasoca-------------------------------------------------
!        if ( noncol ) then
!              call wd_k_points_noncl()
!           else
!              call wd_k_points()
!           endif
! !--------------------------------------------------------------
         e_mpi = 0.d0
         n_mpi = 0
         if(map_k(ik) == myrank_k) then
***************
*** 1934,1942 ****
            end if
            do ie = 1, neg
  !!$             if(map_ek(ie,ik) == mype) n_mpi(ie) = neordr(ie,ik)
!              n_mpi(ie) = neordr(ie,ik)
!              if(map_e(ie) /= myrank_e) cycle
!              e_mpi(ie) = eko_l(map_z(ie),ik)
            end do
         end if
         if(npes>=2) call mpi_allreduce(MPI_IN_PLACE,e_mpi,neg,mpi_double_precision,mpi_sum,mpi_comm_group,ierr)
--- 1940,1958 ----
            end if
            do ie = 1, neg
  !!$             if(map_ek(ie,ik) == mype) n_mpi(ie) = neordr(ie,ik)
! !-----------------------------Modified by T.A.Ariasoca-------------------------------------------------
!              if ( noncol ) then
! 		        jj = ik-1
! 		        ii = ik-mod(jj,2)
!                 n_mpi(ie) = neordr(ie,ii)
!                 if(map_e(ie) /= myrank_e) cycle
!                 e_mpi(ie) = eko_l(map_z(ie),ii)
!              else
!                 n_mpi(ie) = neordr(ie,ik)
!                 if(map_e(ie) /= myrank_e) cycle
!                 e_mpi(ie) = eko_l(map_z(ie),ik)
!              endif
! !--------------------------------------------------------------
            end do
         end if
         if(npes>=2) call mpi_allreduce(MPI_IN_PLACE,e_mpi,neg,mpi_double_precision,mpi_sum,mpi_comm_group,ierr)
***************
*** 2021,2032 ****
           end if
        end if
      end subroutine wd_k_points
    end subroutine m_ESIO_wd_Psicoef
  
    subroutine m_ESIO_wd_BandSymInput(ipri,nfout,nf)
      integer, intent(in) :: ipri,nfout, nf
! 
!     integer :: ik, ie, ri,  nel, ig, ib,ib1,ib2,ibt,ibsize
      integer, allocatable, dimension(:) :: n_mpi
      real(DP),allocatable, dimension(:) :: e_mpi
  
--- 2037,2055 ----
           end if
        end if
      end subroutine wd_k_points
+ !-----------------------------Modified by T.A.Ariasoca-------------------------------------------------
+     subroutine wd_k_points_noncl()
+       if(mype == 0) then
+          write(nf,'(" ik = ",i6,"    ( ",3f14.6," )")') ik,(vkxyz(ik,1:3,BUCS))
+       end if
+     end subroutine wd_k_points_noncl
+ !------------------------------------------------------------------------------
    end subroutine m_ESIO_wd_Psicoef
  
    subroutine m_ESIO_wd_BandSymInput(ipri,nfout,nf)
      integer, intent(in) :: ipri,nfout, nf
!     integer :: ii, jj
!     integer :: ik, ie, ri,  nel, ig, ib,ib1,ib2,ibt,ibsize,asize,asize2
      integer, allocatable, dimension(:) :: n_mpi
      real(DP),allocatable, dimension(:) :: e_mpi
  
***************
*** 2035,2060 ****
  
      allocate(e_mpi(neg)); e_mpi = 0.d0
      allocate(n_mpi(neg)); n_mpi = 0
- 
      if(mype == 0) then
         write(nf,'("##PSIINPSTART")')
         write(nf,'("#MAGNETIC_STATE")')
!        if(nspin == 2 .and. af == 0) then
!           write(nf,'("2  0 ! nspin = 2, af = 0 (FERRO)"/"#")')
         else if(nspin == 2 .and. af == 1) then
!           write(nf,'("2  1 ! nspin = 2, af = 1 (ANTIFERRO)"/"#")')
         else if(nspin == 1) then
!           write(nf,'("1  0 ! nspin = 1, af = 0 (PARAMAGNETIC)"/"#")')
         else
            write(nf,'(2i5," ! magnetic state = unknown"/"#")') nspin, af
         end if
      end if
  
      KPOINT: do ik = 1, kv3, af+1
  
         nel = min(Nw_Psicoef,iba(ik))
         allocate(wf_l(nel,kimg)); wf_l = 0.d0
!        call wd_k_points()
         ! --- Eigen Energies --->
         e_mpi = 0.d0
         n_mpi = 0
--- 2058,2093 ----
  
      allocate(e_mpi(neg)); e_mpi = 0.d0
      allocate(n_mpi(neg)); n_mpi = 0
      if(mype == 0) then
         write(nf,'("##PSIINPSTART")')
         write(nf,'("#MAGNETIC_STATE")')
! 
! !-----------------------------Modified by T.A.Ariasoca-------------------------------------------------
!        if( noncol ) then
!           write(nf,'("2  1  2 ! nspin = 2, af = 1, ndim_spinor = 2 (NONCOLLINEAR)"/"#")')
! 	      else if(nspin == 2 .and. af == 0) then
!           write(nf,'("2  0  1 ! nspin = 2, af = 0, ndim_spinor = 1 (FERRO)"/"#")')
         else if(nspin == 2 .and. af == 1) then
!           write(nf,'("2  1  1 ! nspin = 2, af = 1, ndim_spinor = 1 (ANTIFERRO)"/"#")')
         else if(nspin == 1) then
!           write(nf,'("1  0  1 ! nspin = 1, af = 0, ndim_spinor = 1 (PARAMAGNETIC)"/"#")')
         else
            write(nf,'(2i5," ! magnetic state = unknown"/"#")') nspin, af
         end if
+ !-----------------------------------------------------------------------------------
      end if
  
      KPOINT: do ik = 1, kv3, af+1
  
         nel = min(Nw_Psicoef,iba(ik))
         allocate(wf_l(nel,kimg)); wf_l = 0.d0
! !-----------------------------Modified by T.A.Ariasoca-------------------------------------------------
!        if ( noncol ) then
!              call wd_k_points_noncl()
!           else
!              call wd_k_points()
!           endif
! !----------------------------------------------
         ! --- Eigen Energies --->
         e_mpi = 0.d0
         n_mpi = 0
***************
*** 2065,2073 ****
            end if
            do ie = 1, neg
  !!$             if(map_ek(ie,ik) == mype) n_mpi(ie) = neordr(ie,ik)
!              n_mpi(ie) = neordr(ie,ik)
!              if(map_e(ie) /= myrank_e) cycle
!              e_mpi(ie) = eko_l(map_z(ie),ik)
            end do
         end if
         if(npes>=2) call mpi_allreduce(MPI_IN_PLACE,e_mpi,neg,mpi_double_precision,mpi_sum,mpi_comm_group,ierr)
--- 2098,2116 ----
            end if
            do ie = 1, neg
  !!$             if(map_ek(ie,ik) == mype) n_mpi(ie) = neordr(ie,ik)
! !-----------------------------Modified by T.A.Ariasoca-------------------------------------------------
!              if ( noncol ) then
! 		        jj = ik-1
! 		        ii = ik-mod(jj,2)
!                 n_mpi(ie) = neordr(ie,ii)
!                 if(map_e(ie) /= myrank_e) cycle
!                 e_mpi(ie) = eko_l(map_z(ie),ii)
!              else
!                 n_mpi(ie) = neordr(ie,ik)
!                 if(map_e(ie) /= myrank_e) cycle
!                 e_mpi(ie) = eko_l(map_z(ie),ik)
!              endif
! !-----------------------------------------------
            end do
         end if
         if(npes>=2) call mpi_allreduce(MPI_IN_PLACE,e_mpi,neg,mpi_double_precision,mpi_sum,mpi_comm_group,ierr)
***************
*** 2160,2165 ****
--- 2203,2214 ----
           write(nf,'("#KPOINT"/,i06,3f14.6,2x,a4,/"#")') ik,(vkxyz(ik,1:3,BUCS)),spinstate(ip)
        end if
      end subroutine wd_k_points
+ !-----------------------------Modified by T.A.Ariasoca-------------------------------------------------
+     subroutine wd_k_points_noncl()
+       write(nf,'("#KPOINT"/,i06,3f14.6,2x,a4,/"#")') ik,(vkxyz(ik,1:3,BUCS))
+     end subroutine wd_k_points_noncl
+ !------------------------------------------------------------------------------
+ 
    end subroutine m_ESIO_wd_BandSymInput
  
    subroutine m_ESIO_wd_WFs_and_EVs_ek(nfout,nf)
diff -rcN phase0_2019.01/src_phase/m_ES_dos.F90 phase0_2019.02/src_phase/m_ES_dos.F90
*** phase0_2019.01/src_phase/m_ES_dos.F90	2019-03-19 15:13:12.000000000 +0900
--- phase0_2019.02/src_phase/m_ES_dos.F90	2019-07-26 09:11:10.216144072 +0900
***************
*** 883,889 ****
      
      if(ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if(sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      do ispin = 1, nspin
--- 883,889 ----
      
      if(ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if(iwsc==TOTAL .and. sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      do ispin = 1, nspin
***************
*** 959,965 ****
         do id = 1, nEWindows-1
            sumdos(id+1,ispin) = sumdos(id,ispin) + dos(id+1,ispin)*DeltaE
         end do
!        if(sw_pdos == ON) then
            do iorb = 1,nlmta_phi
               sumpdos(1,iorb,ispin) = pdos(1,iorb,ispin)*DeltaE
               do id = 1, nEWindows-1
--- 959,965 ----
         do id = 1, nEWindows-1
            sumdos(id+1,ispin) = sumdos(id,ispin) + dos(id+1,ispin)*DeltaE
         end do
!        if(iwsc == TOTAL .and. sw_pdos == ON) then
            do iorb = 1,nlmta_phi
               sumpdos(1,iorb,ispin) = pdos(1,iorb,ispin)*DeltaE
               do id = 1, nEWindows-1
***************
*** 999,1005 ****
      
      if (ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if (sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      
--- 999,1005 ----
      
      if (ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if (iwsc == TOTAL .and. sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      
***************
*** 1088,1094 ****
         sumdos(id+1,:) = sumdos(id,:) + dos(id+1,:)*DeltaE
      end do
  
!     if(sw_pdos == ON) then
         do iorb = 1,num_orbitals
            sumpdos(1,iorb,:) = pdos(1,iorb,:)*DeltaE
            do id = 1, nEWindows-1
--- 1088,1094 ----
         sumdos(id+1,:) = sumdos(id,:) + dos(id+1,:)*DeltaE
      end do
  
!     if(iwsc == TOTAL .and. sw_pdos == ON) then
         do iorb = 1,num_orbitals
            sumpdos(1,iorb,:) = pdos(1,iorb,:)*DeltaE
            do id = 1, nEWindows-1
***************
*** 1435,1441 ****
  
      if(ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if(sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      do ispin = 1, nspin
--- 1435,1441 ----
  
      if(ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if(iwsc==TOTAL .and. sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      do ispin = 1, nspin
***************
*** 1506,1512 ****
         do id = 1, nEWindows-1
            sumdos(id+1,ispin) = sumdos(id,ispin) + dos(id+1,ispin)*DeltaE
         end do
!        if(sw_pdos == ON) then
            do iorb = 1,nlmta_phi
               sumpdos(1,iorb,ispin) = pdos(1,iorb,ispin)*DeltaE
               do id = 1, nEWindows-1
--- 1506,1512 ----
         do id = 1, nEWindows-1
            sumdos(id+1,ispin) = sumdos(id,ispin) + dos(id+1,ispin)*DeltaE
         end do
!        if(iwsc == TOTAL .and. sw_pdos == ON) then
            do iorb = 1,nlmta_phi
               sumpdos(1,iorb,ispin) = pdos(1,iorb,ispin)*DeltaE
               do id = 1, nEWindows-1
***************
*** 2238,2244 ****
  !!$                write(nfout,'(" !dos: ",8f9.5)') (eeig2(ip2,ib,ispin),ib=1,neg)
               end do
            end do
!           if(sw_pdos == ON) then
            do ispin=1,nspin
               write(nfout,'(" !dos: ispin = ",i5)') ispin
               do ip2 = 1, np2 
--- 2238,2244 ----
  !!$                write(nfout,'(" !dos: ",8f9.5)') (eeig2(ip2,ib,ispin),ib=1,neg)
               end do
            end do
!           if(icomponent == TOTAL .and. sw_pdos == ON) then
            do ispin=1,nspin
               write(nfout,'(" !dos: ispin = ",i5)') ispin
               do ip2 = 1, np2 
***************
*** 2705,2711 ****
            end do
         end do
  
!        if (sw_pdos == ON) then
  
            do ip2 = 1, np2 
               ik = ndim_spinor *(ip2-1) + 1
--- 2705,2711 ----
            end do
         end do
  
!        if (icomponent == TOTAL .and. sw_pdos == ON) then
  
            do ip2 = 1, np2 
               ik = ndim_spinor *(ip2-1) + 1
diff -rcN phase0_2019.01/src_phase/m_Ionic_System.F90 phase0_2019.02/src_phase/m_Ionic_System.F90
*** phase0_2019.01/src_phase/m_Ionic_System.F90	2019-04-09 14:47:44.000000000 +0900
--- phase0_2019.02/src_phase/m_Ionic_System.F90	2019-08-30 13:39:15.165052056 +0900
***************
*** 12162,12168 ****
         integer :: i,j
         real(kind=DP) :: w,z,c6,n1,n2,dtmp,dtmpr,dwr,dzr
         real(kind=DP), dimension(3) :: ddtmp,dw,dz
!        real(kind=DP) :: eps=1.d-10
         z = 0.d0;w=0.d0;dw=0.d0;dz=0.d0;dwr=0.d0;dzr=0.d0
         do i=1,dftd3par%maxnc(ielem)
            do j=1,dftd3par%maxnc(jelem)
--- 12162,12169 ----
         integer :: i,j
         real(kind=DP) :: w,z,c6,n1,n2,dtmp,dtmpr,dwr,dzr
         real(kind=DP), dimension(3) :: ddtmp,dw,dz
! !       real(kind=DP) :: eps=1.d-10
!        real(kind=DP) :: eps=0.d0
         z = 0.d0;w=0.d0;dw=0.d0;dz=0.d0;dwr=0.d0;dzr=0.d0
         do i=1,dftd3par%maxnc(ielem)
            do j=1,dftd3par%maxnc(jelem)
diff -rcN phase0_2019.01/src_phase/m_Phonon.F90 phase0_2019.02/src_phase/m_Phonon.F90
*** phase0_2019.01/src_phase/m_Phonon.F90	2019-04-02 14:41:39.000000000 +0900
--- phase0_2019.02/src_phase/m_Phonon.F90	2019-07-26 09:26:46.769336142 +0900
***************
*** 2091,2100 ****
      if(phonon_method==PHONON_DOS .and. way_of_smearing==FERMI_DIRAC) free_e=0.d0
      do iq=1,nqvec
         if(phonon_method == PHONON_DOS) then
!           write(nfmode,'(1x,"iq=",i5," q=(",f10.5,",",f10.5,",",f10.5,") (",f10.5,",",f10.5,",",f10.5,")",1x,f15.10)') &
                 & iq,qvin(iq,1:3),qvec(iq,1:3),wght(iq)
         else if(phonon_method /= PHONON_GAMMA) then
!           write(nfmode,'(1x,"iq=",i5," q=(",f10.5,",",f10.5,",",f10.5,") (",f10.5,",",f10.5,",",f10.5,")")') &
                 & iq,qvin(iq,1:3),qvec(iq,1:3)
         end if
         do i=1,nmodes
--- 2091,2100 ----
      if(phonon_method==PHONON_DOS .and. way_of_smearing==FERMI_DIRAC) free_e=0.d0
      do iq=1,nqvec
         if(phonon_method == PHONON_DOS) then
!           write(nfmode,'(1x,"iq=",i5," q=(",f15.10,",",f15.10,",",f15.10,") (",f15.10,",",f15.10,",",f15.10,")",1x,f15.10)') &
                 & iq,qvin(iq,1:3),qvec(iq,1:3),wght(iq)
         else if(phonon_method /= PHONON_GAMMA) then
!           write(nfmode,'(1x,"iq=",i5," q=(",f15.10,",",f15.10,",",f15.10,") (",f15.10,",",f15.10,",",f15.10,")")') &
                 & iq,qvin(iq,1:3),qvec(iq,1:3)
         end if
         do i=1,nmodes
***************
*** 2288,2294 ****
        if(mype /= 0) return
        call m_Files_open_nfphdos()
  
!       write(nfphdos,'("# Index Omega(mHa) Omega(eV) Omega(cm-1) DOS(States/Ha) DOS(States/eV) DOS(States/cm-1) IntDOS(States)")') 
        do ie = 0, newin
           eev = e(ie)*ha2ev/1.0d3  ! Hartree
           ecminv = eev*ev2cminv
--- 2288,2294 ----
        if(mype /= 0) return
        call m_Files_open_nfphdos()
  
!       write(nfphdos,'("# Index Omega(mHa) Omega(eV) Omega(cm-1) DOS(States/mHa) DOS(States/eV) DOS(States/cm-1) IntDOS(States)")') 
        do ie = 0, newin
           eev = e(ie)*ha2ev/1.0d3  ! Hartree
           ecminv = eev*ev2cminv
diff -rcN phase0_2019.01/src_phase/m_PseudoPotential.F90 phase0_2019.02/src_phase/m_PseudoPotential.F90
*** phase0_2019.01/src_phase/m_PseudoPotential.F90	2019-03-19 15:13:12.000000000 +0900
--- phase0_2019.02/src_phase/m_PseudoPotential.F90	2019-07-26 08:12:24.610088449 +0900
***************
*** 10563,10568 ****
--- 10563,10569 ----
  !         do iopr=1,nopr
            do iopr=1,nopr+af !ASMS
               ja=napt(ia,iopr)
+              if(ja>natm) ja=ja-natm
               nrorb(ilmta,iopr) = nylm(isph,iopr)
               do mm=1,nylm(isph,iopr)
  ! debug
diff -rcN phase0_2019.01/src_phase_3d/Initial_Electronic_Structure.F90 phase0_2019.02/src_phase_3d/Initial_Electronic_Structure.F90
*** phase0_2019.01/src_phase_3d/Initial_Electronic_Structure.F90	2019-03-19 15:15:39.000000000 +0900
--- phase0_2019.02/src_phase_3d/Initial_Electronic_Structure.F90	2019-07-26 09:09:55.094188535 +0900
***************
*** 40,46 ****
         &                         , by_pseudo_atomic_orbitals &
         &                         , Valence_plus_PC_Charge, VXC_AND_EXC &
         &                         , ONE_BY_ONE, ALL_AT_ONCE, DP, YES &
!        &                         , Partial_Core_Charge
    use m_IterationNumbers,   only : iteration, nk_in_the_process
    use m_Files,              only : nfout,nfchgt,nfzaj,nfcntn_bin,nfeng,nfefermi &
         &                         , F_ZAJ_in_partitioned, F_CHGT_in_partitioned &
--- 40,46 ----
         &                         , by_pseudo_atomic_orbitals &
         &                         , Valence_plus_PC_Charge, VXC_AND_EXC &
         &                         , ONE_BY_ONE, ALL_AT_ONCE, DP, YES &
!        &                         , Partial_Core_Charge, PT_CONTROL
    use m_IterationNumbers,   only : iteration, nk_in_the_process
    use m_Files,              only : nfout,nfchgt,nfzaj,nfcntn_bin,nfeng,nfefermi &
         &                         , F_ZAJ_in_partitioned, F_CHGT_in_partitioned &
***************
*** 67,73 ****
         &                         , sw_hybrid_functional, sw_screening_correction &
         &                         , sw_external_potential, sw_fef, initial_occmat,kimg &
         &                         , sw_berry_phase, sw_rsb, sw_eval_epc_on_fftmesh &
!        &                         , sw_initial_es
    use m_Kpoints,            only : kv3
    use m_PlaneWaveBasisSet,  only : kg1, m_pwBS_alloc_ylm_l,kgp,ngabc
    use m_Total_Energy,      only  : m_TE_set_etotal_old
--- 67,73 ----
         &                         , sw_hybrid_functional, sw_screening_correction &
         &                         , sw_external_potential, sw_fef, initial_occmat,kimg &
         &                         , sw_berry_phase, sw_rsb, sw_eval_epc_on_fftmesh &
!        &                         , sw_initial_es, imdalg
    use m_Kpoints,            only : kv3
    use m_PlaneWaveBasisSet,  only : kg1, m_pwBS_alloc_ylm_l,kgp,ngabc
    use m_Total_Energy,      only  : m_TE_set_etotal_old
***************
*** 363,369 ****
       end if
    end if
  
- 
    if(sw_phonon == ON .and. sw_calc_force == OFF) return
  
  !!$  call m_PP_gfqwei(nfout)  ! -> modnrm, fqwei, nlmta1, nlmta2
--- 363,368 ----
***************
*** 413,420 ****
  !                    call mpi_barrier(mpi_comm_group, ierr)
                      call timer_sta(30)
  #endif
! !        call m_CD_initial_CD_by_Gauss_func(nfout)   ! (intchg) -> chgq_l
!         call m_CD_initial_CD_by_Gauss_kt(nfout)   ! (intchg) -> chgq_l
  #ifdef FJ_TIMER
                      call timer_end(30)
  #endif
--- 412,422 ----
  !                    call mpi_barrier(mpi_comm_group, ierr)
                      call timer_sta(30)
  #endif
!         if (iteration_unit_cell > 1 .and. sw_read_nfchgt_prev_cell == ON ) then
!            call read_charge_density(condition=-4)
!         else
!            call m_CD_initial_CD_by_Gauss_kt(nfout)   ! (intchg) -> chgq_l
!         end if
  #ifdef FJ_TIMER
                      call timer_end(30)
  #endif
***************
*** 422,428 ****
       end if
  
       !---- set wave functions ----
!      if ( iteration_unit_cell > 1 .and. (sw_read_nfzaj_prev_cell == ON .or. is_charge_density_read)) then
          call read_zaj( condition =-4 )
       else if(intzaj == by_random_numbers) then
  #ifdef FJ_TIMER
--- 424,431 ----
       end if
  
       !---- set wave functions ----
!      if ( iteration_unit_cell > 1 .and. imdalg.ne.PT_CONTROL&
!           & .and. (sw_read_nfzaj_prev_cell == ON .or. is_charge_density_read)) then
          call read_zaj( condition =-4 )
       else if(intzaj == by_random_numbers) then
  #ifdef FJ_TIMER
***************
*** 703,716 ****
    end subroutine read_occ_mat
  
    subroutine read_zaj( condition )
      integer, intent(in) :: condition
  
      call m_Files_open_nfzaj()
  
      if ( condition == 1 ) then
!           call m_ESIO_rd_WFs(nfout,nfzaj,F_ZAJ_in_partitioned)
  
      else if ( condition == -4 ) then       ! coordinate-continuation
      endif
  
    end subroutine read_zaj
--- 706,721 ----
    end subroutine read_occ_mat
  
    subroutine read_zaj( condition )
+     use m_ES_IO,                   only : m_ESIO_import_WFs_prev_cell
      integer, intent(in) :: condition
  
      call m_Files_open_nfzaj()
  
      if ( condition == 1 ) then
!        call m_ESIO_rd_WFs(nfout,nfzaj,F_ZAJ_in_partitioned)
  
      else if ( condition == -4 ) then       ! coordinate-continuation
+        call m_ESIO_import_WFs_prev_cell(nfout,nfzaj,F_ZAJ_in_partitioned)
      endif
  
    end subroutine read_zaj
diff -rcN phase0_2019.01/src_phase_3d/Preparation.F90 phase0_2019.02/src_phase_3d/Preparation.F90
*** phase0_2019.01/src_phase_3d/Preparation.F90	2019-03-19 15:15:39.000000000 +0900
--- phase0_2019.02/src_phase_3d/Preparation.F90	2019-07-26 09:09:55.090188363 +0900
***************
*** 57,62 ****
--- 57,63 ----
         &                       ,P_CONTROL, PT_CONTROL &
         &                       ,DRIVER_URAMP,DRIVER_SC_DFT
    use m_Parallelization,   only:mpi_comm_group &
+        &                       ,m_Parallel_store_prev_np_g1k &
         &                       ,m_Parallel_init_mpi_kngp_3D &
         &                       ,m_Parallel_init_mpi_kngp_B_3D &
  !F       &                       ,m_Parallel_init_mpi_kngp &
***************
*** 146,152 ****
         &                       ,natm,ityp,iatomn,ntyp &
         &                       ,m_IS_reset_extrpl_status
    use m_PseudoPotential,   only:m_PP_input_xctype,ival,m_PP_renew_etot1
!   use m_Kpoints,           only:way_ksample &
         &                       ,m_Kp_gnrt_or_rd_k_points &
         &                       ,m_Kp_alloc_kpoints &
         &                       ,m_Kp_cr_kpoints_table &
--- 147,153 ----
         &                       ,natm,ityp,iatomn,ntyp &
         &                       ,m_IS_reset_extrpl_status
    use m_PseudoPotential,   only:m_PP_input_xctype,ival,m_PP_renew_etot1
!   use m_Kpoints,           only:way_ksample,kv3 &
         &                       ,m_Kp_gnrt_or_rd_k_points &
         &                       ,m_Kp_alloc_kpoints &
         &                       ,m_Kp_cr_kpoints_table &
***************
*** 274,279 ****
--- 275,281 ----
          call m_ESrmm_dealloc_r_norm_flag()
          call m_ESsubmat_dealloc()
          call m_pwBS_store_prev_kg1_kgp()
+         call m_Parallel_store_prev_np_g1k(kv3)
          call m_FFT_reset_firstcall()
       endif
       call m_IS_gdiis_reset()
diff -rcN phase0_2019.01/src_phase_3d/m_ES_IO.F90 phase0_2019.02/src_phase_3d/m_ES_IO.F90
*** phase0_2019.01/src_phase_3d/m_ES_IO.F90	2019-03-19 15:15:39.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_ES_IO.F90	2019-07-26 09:09:55.084188167 +0900
***************
*** 77,83 ****
         &                           , mype,ierr,map_k, map_ek,ista_e,iend_e,istep_e,map_z, np_e &
         &                           , ista_k,iend_k,myrank_e,myrank_k,map_e,nrank_e &
         &                           , ista_kngp,iend_kngp, nrank_k  &
!        &                           , ista_g1k,iend_g1k, np_g1k , myrank_g, nrank_g
    use m_IterationNumbers,     only : nk_in_the_process, nk_converged, nkgroup &
         &                           , first_kpoint_in_this_job, iteration_ionic, iteration_electronic
    use m_FFT,                  only : fft_box_size_WF,nfft
--- 77,83 ----
         &                           , mype,ierr,map_k, map_ek,ista_e,iend_e,istep_e,map_z, np_e &
         &                           , ista_k,iend_k,myrank_e,myrank_k,map_e,nrank_e &
         &                           , ista_kngp,iend_kngp, nrank_k  &
!        &                           , ista_g1k,iend_g1k, np_g1k , myrank_g, nrank_g, np_g1k_prev
    use m_IterationNumbers,     only : nk_in_the_process, nk_converged, nkgroup &
         &                           , first_kpoint_in_this_job, iteration_ionic, iteration_electronic
    use m_FFT,                  only : fft_box_size_WF,nfft
***************
*** 1134,1140 ****
--- 1134,1254 ----
                                                    __TIMER_SUB_STOP(1372)
    end subroutine m_ESIO_rd_WFs
  
+ ! ==== TY 2019/06/25 revised the same subroutine in the 2D version m_ES_IO.f90
+ ! ==== EXP_CELLOPT ==== 2015/09/24 
+   subroutine m_ESIO_import_WFs_prev_cell(nfout,nfzaj, F_ZAJ_partitioned)
+     integer, intent(in) :: nfout, nfzaj
+     logical, intent(in) :: F_ZAJ_partitioned
+     integer    :: ik,ib,ri, i, j
+     integer    :: id_sname = -1
+     integer    :: ierror
+ 
+     call tstatc0_begin('m_ESIO_import_WFs_prev_cell ',id_sname)
  
+     if(precision_WFfile==SP) then
+        if(ipri >= 1) write(nfout,*) ' !D Reading zaj (single_precision) in m_ESIO_import_WFs_prev_cell'
+     else
+        if(ipri >= 1) write(nfout,*) ' !D Reading zaj (double_precision) in m_ESIO_import_WFs_prev_cell'
+     end if
+     rewind nfzaj
+     if(F_ZAJ_partitioned) then
+        if(precision_WFfile==SP) then
+           allocate(wf_l(maxval(np_g1k_prev),kimg)); wf_l = 0.d0
+        else
+           allocate(wfdp_l(maxval(np_g1k_prev),kimg)); wfdp_l = 0.d0
+        end if
+ !!$
+ !!$    zaj_l = 0.0d0
+        do ik = ista_k, iend_k, af+1        ! MPI
+           do ib = 1, np_e
+ !!$          do ib = ista_e, iend_e, istep_e  ! MPI
+              if(ista_e+ib-1 >neg_previous) cycle
+ !!$             if(ib > neg_previous) cycle
+              if(precision_WFfile==SP) then
+                 read(nfzaj) wf_l
+ 
+                 if(kimg == 1) then
+                    do i = 1, min(np_g1k(ik),np_g1k_prev(ik))
+                       zaj_l(i,ib,ik,1) = wf_l(i,1)
+                    end do
+                 else if(kimg==2) then
+                    do i = 1, min(np_g1k(ik),np_g1k_prev(ik))
+                       zaj_l(i,ib,ik,1) = wf_l(i,1)
+                       zaj_l(i,ib,ik,2) = wf_l(i,2)
+                    end do
+                 end if
+              else if(precision_WFfile==DP) then
+                 read(nfzaj) wfdp_l
+ 
+                 if(kimg == 1) then
+                    do i = 1, min(np_g1k(ik),np_g1k_prev(ik))
+                       zaj_l(i,ib,ik,1) = wfdp_l(i,1)
+                    end do
+                 else if(kimg==2) then
+                    do i = 1, min(np_g1k(ik),np_g1k_prev(ik))
+                       zaj_l(i,ib,ik,1) = wfdp_l(i,1)
+                       zaj_l(i,ib,ik,2) = wfdp_l(i,2)
+                    end do
+                 end if
+ 
+              end if
+           end do
+ 
+        end do
+        if(precision_WFfile==SP) then
+           deallocate(wf_l)
+        else if(precision_WFfile==DP) then
+           deallocate(wfdp_l)
+        end if
+     else
+        zaj_l = 0.d0
+        if(precision_WFfile==SP) then
+           allocate(wf_l(kg1_prev,kimg));  wf_l = 0
+        else
+           allocate(wfdp_l(kg1_prev,kimg));  wfdp_l = 0
+        end if
+        do ik = 1, kv3, af+1
+           do ib = 1, neg_previous
+ ! -----------------
+            if(precision_WFfile==SP) then
+              if(mype == 0) read(nfzaj, end = 9999, err = 9999) wf_l  ! MPI
+              call mpi_bcast(wf_l,kg1_prev*kimg,mpi_real,0,mpi_comm_group,ierr)
+ 
+              if((map_k(ik) == myrank_k) .and. (map_e(ib) == myrank_e) )then         ! MPI
+                 do ri = 1, kimg
+                   do j = min(ista_g1k(ik),kg1_prev),min(iend_g1k(ik),kg1_prev)
+                      zaj_l(j-ista_g1k(ik)+1,map_z(ib),ik,ri) = wf_l(j,ri)  ! MPI
+                   end do
+                 end do
+              end if
+ ! -----------------
+            else if(precision_WFfile==DP) then
+              if(mype == 0) read(nfzaj, end = 9999, err = 9999) wfdp_l  ! MPI
+              call mpi_bcast(wfdp_l,kg1_prev*kimg,mpi_double_precision,0,mpi_comm_group,ierr)
+              if((map_k(ik) == myrank_k) .and. (map_e(ib) == myrank_e) )then         ! MPI
+                 do ri = 1, kimg
+                   do j = min(ista_g1k(ik),kg1_prev),min(iend_g1k(ik),kg1_prev)
+                      zaj_l(j-ista_g1k(ik)+1,map_z(ib),ik,ri) = wfdp_l(j,ri)  ! MPI
+                   end do
+                 end do
+              end if
+            endif
+ ! -----------------
+           end do
+        end do
+        if(precision_WFfile==SP) then
+           deallocate(wf_l)
+        else if(precision_WFfile==DP) then
+           deallocate(wfdp_l)
+        end if
+     end if
+     return
+ 9999 continue
+     ierror = EOF_REACHED
+     call phase_error_wo_filename(ierror, nfout, nfzaj, __LINE__, __FILE__)
+   end subroutine m_ESIO_import_WFs_prev_cell
+ ! ===================== 2015/09/24
+ ! ==== TY 2019/06/25
  
    subroutine m_ESIO_wd_WFs_standardout(nfout,ipriwf)
      integer, intent(in) :: nfout,ipriwf
diff -rcN phase0_2019.01/src_phase_3d/m_ES_WF_by_submat.F90 phase0_2019.02/src_phase_3d/m_ES_WF_by_submat.F90
*** phase0_2019.01/src_phase_3d/m_ES_WF_by_submat.F90	2019-03-19 15:15:38.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_ES_WF_by_submat.F90	2019-07-26 09:09:55.075187763 +0900
***************
*** 457,462 ****
--- 457,466 ----
  #endif
  
      integer :: isrsize, num
+ !**** TYamasaki 2019/06/10
+     integer :: nb,kb,jb,ib
+ !**** <===
+ 
      integer, dimension(0:nrank_g-1) ::req_r,req_s
      integer, dimension(MPI_STATUS_SIZE,0:nrank_g-1)::sta_r, sta_s
       real(kind=DP), allocatable, dimension(:,:) :: sendbuf, recvbuf
***************
*** 518,535 ****
      call mpi_allreduce(MPI_IN_PLACE,eko_d,neg,mpi_double_precision,mpi_sum,mpi_kg_world,ierr)
                                                    __TIMER_COMM_STOP(914)
      eko1 = sum(eko_d(1:neg))
!     allocate(wk1(maxval(np_g1k),neg,kimg))
! ! === DEBUG by tkato 2011/10/24 ================================================
!     wk1 = 0.0d0
!     do i = ista_e, iend_e
!        wk1(:,i,:) = zaj_l(:,i-ista_e+1,ik,:)
!     enddo
!     call mpi_allreduce(MPI_IN_PLACE,wk1,maxval(np_g1k)*neg*kimg,MPI_DOUBLE_PRECISION,MPI_SUM,mpi_kg_world,ierr)
!     do i = 1, neg
!        zaj_ball(:,i,ik,:) = wk1(:,i,:)
!     enddo
  ! ==============================================================================
!     deallocate(wk1)
  
      if(sw_gep == ON)then
         allocate(zat_l(maxval(np_g1k),np_e,ista_k:iend_k,kimg)); zat_l = zaj_l
--- 522,566 ----
      call mpi_allreduce(MPI_IN_PLACE,eko_d,neg,mpi_double_precision,mpi_sum,mpi_kg_world,ierr)
                                                    __TIMER_COMM_STOP(914)
      eko1 = sum(eko_d(1:neg))
! ! == revised by T.Yamasaki 2019/06/10 ==============================================================================
! !!$    allocate(wk1(maxval(np_g1k),neg,kimg))
!     if(nrank_e == 1) then
!        zaj_ball(:,:,ik,:) = zaj_l(:,:,ik,:)
!     else
! !!$       write(nfout,'(" *** mpi_allreduce size = ",i12," *** ")') maxval(np_g1k)*neg*kimg*8
! !!$       write(nfout,'(" ***       reduced size = ",i12," *** ")') maxval(np_g1k)*mp_e*kimg
! !!$       call flush(nfout)
!        allocate(wk1(maxval(np_g1k),mp_e,kimg),stat=ierr)
!        do nb = 0, nrank_e-1
!           if(nb == myrank_e) then
!              wk1(:,1:np_e,1:kimg) = zaj_l(:,1:np_e,ik,1:kimg)
!              if(np_e+1<mp_e) wk1(:,np_e+1:mp_e,:) = 0.d0
!           end if
!           call mpi_bcast(wk1,maxval(np_g1k)*mp_e*kimg,mpi_double_precision,nb,mpi_kg_world,ierr)
!           do kb = 1, kimg
!              do jb = nis_e(nb),nie_e(nb)
!                 do ib = 1, np_g1k(ik)
!                    zaj_ball(ib,jb,ik,kb) = wk1(ib,jb-nis_e(nb)+1,kb)
!                 end do
!              end do
!           end do
!        end do
!        deallocate(wk1)
!     end if
! !!$! === DEBUG by tkato 2011/10/24 ================================================
! !!$    wk1 = 0.0d0
! !!$    write(nfout,'(" *** mpi_allreduce size = ",i12," *** ")') maxval(np_g1k)*neg*kimg*8
! !!$    do i = ista_e, iend_e
! !!$       wk1(:,i,:) = zaj_l(:,i-ista_e+1,ik,:)
! !!$    enddo
! !!$    call mpi_allreduce(MPI_IN_PLACE,wk1,maxval(np_g1k)*neg*kimg,MPI_DOUBLE_PRECISION,MPI_SUM,mpi_kg_world,ierr)
! !!$    do i = 1, neg
! !!$       zaj_ball(:,i,ik,:) = wk1(:,i,:)
! !!$    enddo
! !!$! ==============================================================================
! !!$    deallocate(wk1)
  ! ==============================================================================
! 
  
      if(sw_gep == ON)then
         allocate(zat_l(maxval(np_g1k),np_e,ista_k:iend_k,kimg)); zat_l = zaj_l
diff -rcN phase0_2019.01/src_phase_3d/m_ES_dos.F90 phase0_2019.02/src_phase_3d/m_ES_dos.F90
*** phase0_2019.01/src_phase_3d/m_ES_dos.F90	2019-03-19 15:15:39.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_ES_dos.F90	2019-07-26 09:09:55.090188363 +0900
***************
*** 892,898 ****
      
      if(ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if(sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      do ispin = 1, nspin
--- 892,898 ----
      
      if(ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if(iwsc == TOTAL .and. sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      do ispin = 1, nspin
***************
*** 968,974 ****
         do id = 1, nEWindows-1
            sumdos(id+1,ispin) = sumdos(id,ispin) + dos(id+1,ispin)*DeltaE
         end do
!        if(sw_pdos == ON) then
            do iorb = 1,nlmta_phi
               sumpdos(1,iorb,ispin) = pdos(1,iorb,ispin)*DeltaE
               do id = 1, nEWindows-1
--- 968,974 ----
         do id = 1, nEWindows-1
            sumdos(id+1,ispin) = sumdos(id,ispin) + dos(id+1,ispin)*DeltaE
         end do
!        if(iwsc == TOTAL .and. sw_pdos == ON) then
            do iorb = 1,nlmta_phi
               sumpdos(1,iorb,ispin) = pdos(1,iorb,ispin)*DeltaE
               do id = 1, nEWindows-1
***************
*** 1008,1014 ****
      
      if (ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if (sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      
--- 1008,1014 ----
      
      if (ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if (iwsc == TOTAL .and. sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      
***************
*** 1097,1103 ****
         sumdos(id+1,:) = sumdos(id,:) + dos(id+1,:)*DeltaE
      end do
  
!     if(sw_pdos == ON) then
         do iorb = 1,num_orbitals
            sumpdos(1,iorb,:) = pdos(1,iorb,:)*DeltaE
            do id = 1, nEWindows-1
--- 1097,1103 ----
         sumdos(id+1,:) = sumdos(id,:) + dos(id+1,:)*DeltaE
      end do
  
!     if(iwsc == TOTAL .and. sw_pdos == ON) then
         do iorb = 1,num_orbitals
            sumpdos(1,iorb,:) = pdos(1,iorb,:)*DeltaE
            do id = 1, nEWindows-1
***************
*** 1444,1450 ****
  
      if(ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if(sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      do ispin = 1, nspin
--- 1444,1450 ----
  
      if(ipridos >= 2) write(nfout,'(" !! Es, DeltaE = ",2d20.10)') Es, DeltaE
      dos = 0.d0; sumdos = 0.d0
!     if(iwsc == TOTAL .and. sw_pdos == ON) then
         pdos=0.d0; sumpdos=0.d0
      end if
      do ispin = 1, nspin
***************
*** 1515,1521 ****
         do id = 1, nEWindows-1
            sumdos(id+1,ispin) = sumdos(id,ispin) + dos(id+1,ispin)*DeltaE
         end do
!        if(sw_pdos == ON) then
            do iorb = 1,nlmta_phi
               sumpdos(1,iorb,ispin) = pdos(1,iorb,ispin)*DeltaE
               do id = 1, nEWindows-1
--- 1515,1521 ----
         do id = 1, nEWindows-1
            sumdos(id+1,ispin) = sumdos(id,ispin) + dos(id+1,ispin)*DeltaE
         end do
!        if(iwsc == TOTAL .and. sw_pdos == ON) then
            do iorb = 1,nlmta_phi
               sumpdos(1,iorb,ispin) = pdos(1,iorb,ispin)*DeltaE
               do id = 1, nEWindows-1
***************
*** 2256,2262 ****
  !!$                write(nfout,'(" !dos: ",8f9.5)') (eeig2(ip2,ib,ispin),ib=1,neg)
               end do
            end do
!           if(sw_pdos == ON) then
            do ispin=1,nspin
               write(nfout,'(" !dos: ispin = ",i5)') ispin
               do ip2 = 1, np2 
--- 2256,2262 ----
  !!$                write(nfout,'(" !dos: ",8f9.5)') (eeig2(ip2,ib,ispin),ib=1,neg)
               end do
            end do
!           if(icomponent == TOTAL .and. sw_pdos == ON) then
            do ispin=1,nspin
               write(nfout,'(" !dos: ispin = ",i5)') ispin
               do ip2 = 1, np2 
***************
*** 2729,2735 ****
            end do
         end do
  
!        if (sw_pdos == ON) then
  
            do ip2 = 1, np2 
               ik = ndim_spinor *(ip2-1) + 1
--- 2729,2735 ----
            end do
         end do
  
!        if (icomponent == TOTAL .and. sw_pdos == ON) then
  
            do ip2 = 1, np2 
               ik = ndim_spinor *(ip2-1) + 1
diff -rcN phase0_2019.01/src_phase_3d/m_ES_occup.F90 phase0_2019.02/src_phase_3d/m_ES_occup.F90
*** phase0_2019.01/src_phase_3d/m_ES_occup.F90	2019-03-19 15:15:38.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_ES_occup.F90	2019-07-26 09:09:55.075187763 +0900
***************
*** 143,149 ****
    character(len("occ_up")),private,parameter :: tag_occ_up = "occ_up"
    character(len("occ_down")),private,parameter :: tag_occ_down = "occ_down"
  
!   integer, private,parameter :: PRINTOUTLEVEL_TOTCH_SPIN = 2
  
    include 'mpif.h'
  contains
--- 143,150 ----
    character(len("occ_up")),private,parameter :: tag_occ_up = "occ_up"
    character(len("occ_down")),private,parameter :: tag_occ_down = "occ_down"
  
! !!$  integer, private,parameter :: PRINTOUTLEVEL_TOTCH_SPIN = 2
!   integer, private,parameter :: PRINTOUTLEVEL_TOTCH_SPIN = 1
  
    include 'mpif.h'
  contains
***************
*** 1548,1556 ****
  !          call get_occup_l_and_tot(efermi)       ! -(contained here) ->(occup_l,tot)
  !
            if ( noncol ) then
!             call get_occup_l_and_tot_noncl_fix(efermi)     
            else
!             call get_occup_l_and_tot(efermi)     
            endif
  ! ======================================================================= 11.0
  
--- 1549,1557 ----
  !          call get_occup_l_and_tot(efermi)       ! -(contained here) ->(occup_l,tot)
  !
            if ( noncol ) then
!             call get_occup_l_and_tot_noncl_fix(efermi)
            else
!             call get_occup_l_and_tot(efermi)
            endif
  ! ======================================================================= 11.0
  
***************
*** 1790,1796 ****
      jcount = 1
      occupied_ch_equals_totch : do
  
!        call get_occup_l_and_tot(1,1,efermi)      ! -(contained here) ->(occup_l,tot)
  !          ~~~~~~~~~~~~~~~~~~~~
         if(jcount == 1 .and. tot < totch) &
              call wd_fermi_error1(nfout,emin,emax,tot,totch) ! -(b_Fermi)
--- 1791,1797 ----
      jcount = 1
      occupied_ch_equals_totch : do
  
!        call get_occup_l_and_tot(1,1,efermi,.true.)          ! -(contained here) ->(occup_l,tot)
  !          ~~~~~~~~~~~~~~~~~~~~
         if(jcount == 1 .and. tot < totch) &
              call wd_fermi_error1(nfout,emin,emax,tot,totch) ! -(b_Fermi)
***************
*** 1836,1842 ****
            jcount = 1
            occupied_ch_equals_totch2 : do
  
!              call get_occup_l_and_tot(nspin,is,efermi_spin(is))     ! -(contained here) ->(occup_l,tot)
               !    ~~~~~~~~~~~~~~~~~~~
               if(jcount == 1 .and. tot < totch_spin(is)) &
                    call wd_fermi_error1(nfout,emin,emax,tot,totch) ! -(b_Fermi)
--- 1837,1843 ----
            jcount = 1
            occupied_ch_equals_totch2 : do
  
!              call get_occup_l_and_tot(nspin,is,efermi_spin(is),.true.)     ! -(contained here) ->(occup_l,tot)
               !    ~~~~~~~~~~~~~~~~~~~
               if(jcount == 1 .and. tot < totch_spin(is)) &
                    call wd_fermi_error1(nfout,emin,emax,tot,totch) ! -(b_Fermi)
***************
*** 1859,1870 ****
                 & is,efermi_spin(is), tot/2.0
         end do
         call get_total_spin0
-        if(iprioccup >= PRINTOUTLEVEL_TOTCH_SPIN) &
-             & call wd_efermi_and_total_spin0_Plus(nfout,total_spin0,totch_spin0,totch_spin)
         efermi = (efermi_spin(1)+efermi_spin(2))*0.5d0
      else if(imag == FERRO .and. sw_fix_total_spin == NO .and. nspin == 2) then
         call get_total_spin0
-        if(iprioccup >= PRINTOUTLEVEL_TOTCH_SPIN) call wd_efermi_and_total_spin0(nfout,total_spin0,totch_spin0)
      end if
  
  !   do ik = ista_k, iend_k                              ! MPI
--- 1860,1869 ----
                 & is,efermi_spin(is), tot/2.0
         end do
         call get_total_spin0
         efermi = (efermi_spin(1)+efermi_spin(2))*0.5d0
+        if(iprioccup>=1) write(nfout,'(a,f10.4)') ' == new efermi = ',efermi
      else if(imag == FERRO .and. sw_fix_total_spin == NO .and. nspin == 2) then
         call get_total_spin0
      end if
  
  !   do ik = ista_k, iend_k                              ! MPI
***************
*** 1891,1902 ****
      subroutine get_total_spin0
                                                    __TIMER_SUB_START(708)
        do is = 1, nspin
!          call get_occup_l_and_tot(nspin,is,efermi)
           totch_spin0(is) = tot
        end do
!       total_spin0 = totch_spin0(1) - totch_spin0(2)
!       if(iprioccup >= 1) then
!          if(sw_fix_total_spin == YES ) then
              call wd_efermi_and_total_spin0_Plus(nfout,total_spin0,totch_spin0,totch_spin)
           else
              call wd_efermi_and_total_spin0(nfout,total_spin0,totch_spin0)
--- 1890,1901 ----
      subroutine get_total_spin0
                                                    __TIMER_SUB_START(708)
        do is = 1, nspin
!          call get_occup_l_and_tot(nspin,is,efermi,.false.)
           totch_spin0(is) = tot
        end do
!       total_spin0 = (totch_spin0(1) - totch_spin0(2))*0.5d0
!       if(iprioccup >= PRINTOUTLEVEL_TOTCH_SPIN) then
!          if(sw_fix_total_spin == YES) then
              call wd_efermi_and_total_spin0_Plus(nfout,total_spin0,totch_spin0,totch_spin)
           else
              call wd_efermi_and_total_spin0(nfout,total_spin0,totch_spin0)
***************
*** 1905,1915 ****
                                                    __TIMER_SUB_STOP(708)
      end subroutine get_total_spin0
  
!     subroutine get_occup_l_and_tot(nspin,is,efermi)
        integer, intent(in) :: nspin, is
        real(kind=DP), intent(in) :: efermi
!       integer       :: k, i
!       real(kind=DP) :: wspin = 1.d0, e, dos, weight
                                                    __TIMER_SUB_START(704)
        tot = 0.d0
                                                    __TIMER_SUB_START(705)
--- 1904,1915 ----
                                                    __TIMER_SUB_STOP(708)
      end subroutine get_total_spin0
  
!     subroutine get_occup_l_and_tot(nspin,is,efermi,update_occ)
        integer, intent(in) :: nspin, is
        real(kind=DP), intent(in) :: efermi
!       logical, intent(in) :: update_occ
!       integer       :: k, i,iupdown
!       real(kind=DP) :: wspin = 1.d0, e, dos, weight, tmp
                                                    __TIMER_SUB_START(704)
        tot = 0.d0
                                                    __TIMER_SUB_START(705)
***************
*** 1920,1927 ****
           do i = 1, neg
              e = eko_mpi(i,k)
              call width2(e,efermi,width,dos,weight)  ! -(b_Fermi)
!             occup_mpi(i,k) = weight*wspin*kv3*qwgt(k)
!             tot = tot + 2*occup_mpi(i,k)
           end do
        end do
                                                    __TIMER_DO_STOP(804)
--- 1920,1938 ----
           do i = 1, neg
              e = eko_mpi(i,k)
              call width2(e,efermi,width,dos,weight)  ! -(b_Fermi)
!             if(sw_manual_occupation==ON.and..not.update_occ)then
!                if(band_index(i)>0)then
!                   iupdown = 1
!                   if(nspin>1.and.mod(k,2)==0) iupdown = 2
!                   weight = occ_ext(band_index(i),iupdown)
!                endif
!             endif
! ! ====================================================================== 12.1
! 
!             tmp = weight*wspin*kv3*qwgt(k)
!             if(update_occ) occup_mpi(i,k) = tmp
!             !tot = tot + 2*occup_mpi(i,k)
!             tot = tot + 2*tmp
           end do
        end do
                                                    __TIMER_DO_STOP(804)
***************
*** 2345,2361 ****
      call tstatc0_end(id_sname)
                                                    __TIMER_SUB_STOP(715)
    contains
!     subroutine get_total_spin0
        do is = 1, nspin
!          call get_occup_l_and_tot(nspin,is,efermi)
           totch_spin0(is) = tot
        end do
        total_spin0 = totch_spin0(1) - totch_spin0(2)
      end subroutine get_total_spin0
  
!     subroutine get_occup_l_and_tot(nspin,is,efermi)
        integer, intent(in) :: nspin, is
        real(kind=DP), intent(in) :: efermi
        integer       :: k, i
        real(kind=DP) :: wspin = 1.d0, e, dos, weight
        tot = 0.d0
--- 2356,2376 ----
      call tstatc0_end(id_sname)
                                                    __TIMER_SUB_STOP(715)
    contains
!     subroutine get_total_spin0(sw_occupation)
!       integer, optional, intent(in) :: sw_occupation
!       integer :: sw_occupation_t 
!       sw_occupation_t = OFF
        do is = 1, nspin
!          call get_occup_l_and_tot(nspin,is,efermi,sw_occupation_t)
           totch_spin0(is) = tot
        end do
        total_spin0 = totch_spin0(1) - totch_spin0(2)
      end subroutine get_total_spin0
  
!     subroutine get_occup_l_and_tot(nspin,is,efermi,sw_occupation_t)
        integer, intent(in) :: nspin, is
        real(kind=DP), intent(in) :: efermi
+       integer,optional,intent(in) :: sw_occupation_t
        integer       :: k, i
        real(kind=DP) :: wspin = 1.d0, e, dos, weight
        tot = 0.d0
***************
*** 2364,2371 ****
           do i = 1, neg
              e = eko_mpi(i,k)
              call coldsmearing(e,efermi,width,dos,weight)
!             occup_mpi(i,k) = weight*wspin*kv3*qwgt(k)
!             tot = tot + 2*occup_mpi(i,k)
           end do
        end do
        if(af == 1) then
--- 2379,2386 ----
           do i = 1, neg
              e = eko_mpi(i,k)
              call coldsmearing(e,efermi,width,dos,weight)
!             tot = tot + 2*weight*wspin*kv3*qwgt(k)
!             if(.not.(present(sw_occupation_t).and.sw_occupation_t == OFF)) occup_mpi(i,k) = weight*wspin*kv3*qwgt(k)
           end do
        end do
        if(af == 1) then
***************
*** 2505,2520 ****
      integer, intent(in) :: nfout
      real(kind=DP), intent(in) :: total_spin0, totch_spin0(2)
      write(nfout,'(" == efermi = ",f10.4, ", total_spin0 = ",f12.6, ", totch_spin0(1:2) = ",2f12.6)') &
!          & efermi, total_spin0, totch_spin0(1:2)
    end subroutine wd_efermi_and_total_spin0
  
    subroutine wd_efermi_and_total_spin0_plus(nfout,total_spin0,totch_spin0,totch_spin)
      integer, intent(in) :: nfout
      real(kind=DP), intent(in) :: total_spin0, totch_spin0(2),totch_spin(2)
!     write(nfout,'(" == efermi_spin(1:2) = ",2f10.4, ", total_spin  = ",f12.6,", totch_spin(1:2) = ",2f12.6)') &
           & efermi_spin(1:2),total_spin, totch_spin(1:2)*0.5d0
      write(nfout,'(" == efermi = ",f10.4, 20x,", total_spin0 = ",f12.6, ", totch_spin0(1:2) = ",2f12.6)') &
!          & efermi, total_spin0, totch_spin0(1:2)
    end subroutine wd_efermi_and_total_spin0_plus
  
  ! =============== KT_add ========================== 13.0E
--- 2520,2537 ----
      integer, intent(in) :: nfout
      real(kind=DP), intent(in) :: total_spin0, totch_spin0(2)
      write(nfout,'(" == efermi = ",f10.4, ", total_spin0 = ",f12.6, ", totch_spin0(1:2) = ",2f12.6)') &
!          & efermi, total_spin0*0.5d0, totch_spin0(1:2)*0.5d0
    end subroutine wd_efermi_and_total_spin0
  
    subroutine wd_efermi_and_total_spin0_plus(nfout,total_spin0,totch_spin0,totch_spin)
      integer, intent(in) :: nfout
      real(kind=DP), intent(in) :: total_spin0, totch_spin0(2),totch_spin(2)
!     write(nfout,'(" == efermi_spin(1:2) = ",2f10.4, ", total_spin  = ",f12.6,", totch_spin(1:2)  = ",2f12.6)') &
           & efermi_spin(1:2),total_spin, totch_spin(1:2)*0.5d0
+ ! === TY revised 2019.06.13 ======
      write(nfout,'(" == efermi = ",f10.4, 20x,", total_spin0 = ",f12.6, ", totch_spin0(1:2) = ",2f12.6)') &
!          & efermi, total_spin0*0.5d0, totch_spin0(1:2)*0.5d0
! ! ================================
    end subroutine wd_efermi_and_total_spin0_plus
  
  ! =============== KT_add ========================== 13.0E
diff -rcN phase0_2019.01/src_phase_3d/m_ES_ortho.F90 phase0_2019.02/src_phase_3d/m_ES_ortho.F90
*** phase0_2019.01/src_phase_3d/m_ES_ortho.F90	2019-03-19 15:15:38.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_ES_ortho.F90	2019-07-26 09:09:55.074187702 +0900
***************
*** 586,592 ****
            if(nb == myrank_e) then
               do kb = 1, kimg
                  do jb = 1, np_e
!                    wk_zaj(:,jb,kb) = zaj_l(ib,jb,ik,kb)
                  end do
               end do
               if(np_e+1 < mp_e) wk_zaj(:,np_e+1:mp_e,:) = 0.d0
--- 586,592 ----
            if(nb == myrank_e) then
               do kb = 1, kimg
                  do jb = 1, np_e
!                    wk_zaj(:,jb,kb) = zaj_l(:,jb,ik,kb)
                  end do
               end do
               if(np_e+1 < mp_e) wk_zaj(:,np_e+1:mp_e,:) = 0.d0
diff -rcN phase0_2019.01/src_phase_3d/m_Ionic_System.F90 phase0_2019.02/src_phase_3d/m_Ionic_System.F90
*** phase0_2019.01/src_phase_3d/m_Ionic_System.F90	2019-03-20 09:38:56.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_Ionic_System.F90	2019-08-30 13:39:45.841261327 +0900
***************
*** 11429,11435 ****
         integer :: i,j
         real(kind=DP) :: w,z,c6,n1,n2,dtmp,dtmpr,dwr,dzr
         real(kind=DP), dimension(3) :: ddtmp,dw,dz
!        real(kind=DP) :: eps=1.d-10
         z = 0.d0;w=0.d0;dw=0.d0;dz=0.d0;dwr=0.d0;dzr=0.d0
         do i=1,dftd3par%maxnc(ielem)
            do j=1,dftd3par%maxnc(jelem)
--- 11429,11436 ----
         integer :: i,j
         real(kind=DP) :: w,z,c6,n1,n2,dtmp,dtmpr,dwr,dzr
         real(kind=DP), dimension(3) :: ddtmp,dw,dz
! !       real(kind=DP) :: eps=1.d-10
!        real(kind=DP) :: eps=0.d0
         z = 0.d0;w=0.d0;dw=0.d0;dz=0.d0;dwr=0.d0;dzr=0.d0
         do i=1,dftd3par%maxnc(ielem)
            do j=1,dftd3par%maxnc(jelem)
diff -rcN phase0_2019.01/src_phase_3d/m_Parallelization.F90 phase0_2019.02/src_phase_3d/m_Parallelization.F90
*** phase0_2019.01/src_phase_3d/m_Parallelization.F90	2019-03-19 15:15:39.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_Parallelization.F90	2019-07-26 09:09:55.086188230 +0900
***************
*** 189,195 ****
  !!$#ifdef TRANSPOSE
    integer                            :: ista_g1, iend_g1, np_g1, mp_g1
    integer, allocatable, dimension(:) :: nis_g1, nie_g1, nel_g1      !d(0:nrank_g1-1)
!   integer, allocatable, dimension(:) :: ista_g1k, iend_g1k, np_g1k, mp_g1k !d(kv3)
    integer, allocatable, dimension(:,:):: nis_g1k, nie_g1k, nel_g1k  !d(0:nrank_g1-1,kv3)
  
    integer                            :: ista_fs, iend_fs, np_fs, mp_fs
--- 189,195 ----
  !!$#ifdef TRANSPOSE
    integer                            :: ista_g1, iend_g1, np_g1, mp_g1
    integer, allocatable, dimension(:) :: nis_g1, nie_g1, nel_g1      !d(0:nrank_g1-1)
!   integer, allocatable, dimension(:) :: ista_g1k, iend_g1k, np_g1k, mp_g1k, np_g1k_prev !d(kv3)
    integer, allocatable, dimension(:,:):: nis_g1k, nie_g1k, nel_g1k  !d(0:nrank_g1-1,kv3)
  
    integer                            :: ista_fs, iend_fs, np_fs, mp_fs
***************
*** 2046,2051 ****
--- 2046,2052 ----
      if(allocated(ista_g1k)) deallocate(ista_g1k)
      if(allocated(iend_g1k)) deallocate(iend_g1k)
      if(allocated(np_g1k)) deallocate(np_g1k)
+     if(allocated(np_g1k_prev)) deallocate(np_g1k_prev)
      if(allocated(mp_g1k)) deallocate(mp_g1k)
      
      if(allocated(nis_fftp)) deallocate(nis_fftp)
***************
*** 4001,4006 ****
--- 4002,4008 ----
      mpi_nval_enabled = .false.
    end subroutine m_Parallel_dealloc_mpi_nval
  
+ 
    subroutine m_Parallel_init_mpi_iba_3D(nfout,ipri,printable,kv3,iba)
      integer, intent(in) :: nfout,ipri, kv3, iba(kv3)
      logical, intent(in) :: printable
***************
*** 9025,9028 ****
--- 9027,9044 ----
    end subroutine m_Parallel_dealloc_mpi_paw_3D
  #endif
  
+ !! TY 2019.06.25 -->
+   subroutine m_Parallel_store_prev_np_g1k(kv3)
+     integer, intent(in) :: kv3
+     call m_Parallel_alloc_np_g1k_prev(kv3)
+     np_g1k_prev = np_g1k
+   end subroutine m_Parallel_store_prev_np_g1k
+   
+   subroutine m_Parallel_alloc_np_g1k_prev(kv3)
+     integer, intent(in) :: kv3
+     if(allocated(np_g1k_prev)) deallocate(np_g1k_prev);  allocate(np_g1k_prev(kv3))
+   end subroutine m_Parallel_alloc_np_g1k_prev
+ !! <--
+ 
+ 
  end module m_Parallelization
diff -rcN phase0_2019.01/src_phase_3d/m_Phonon.F90 phase0_2019.02/src_phase_3d/m_Phonon.F90
*** phase0_2019.01/src_phase_3d/m_Phonon.F90	2019-04-02 14:42:01.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_Phonon.F90	2019-07-26 09:27:15.114256506 +0900
***************
*** 2091,2100 ****
      if(phonon_method==PHONON_DOS .and. way_of_smearing==FERMI_DIRAC) free_e=0.d0
      do iq=1,nqvec
         if(phonon_method == PHONON_DOS) then
!           write(nfmode,'(1x,"iq=",i5," q=(",f10.5,",",f10.5,",",f10.5,") (",f10.5,",",f10.5,",",f10.5,")",1x,f15.10)') &
                 & iq,qvin(iq,1:3),qvec(iq,1:3),wght(iq)
         else if(phonon_method /= PHONON_GAMMA) then
!           write(nfmode,'(1x,"iq=",i5," q=(",f10.5,",",f10.5,",",f10.5,") (",f10.5,",",f10.5,",",f10.5,")")') &
                 & iq,qvin(iq,1:3),qvec(iq,1:3)
         end if
         do i=1,nmodes
--- 2091,2100 ----
      if(phonon_method==PHONON_DOS .and. way_of_smearing==FERMI_DIRAC) free_e=0.d0
      do iq=1,nqvec
         if(phonon_method == PHONON_DOS) then
!           write(nfmode,'(1x,"iq=",i5," q=(",f15.10,",",f15.10,",",f15.10,") (",f15.10,",",f15.10,",",f15.10,")",1x,f15.10)') &
                 & iq,qvin(iq,1:3),qvec(iq,1:3),wght(iq)
         else if(phonon_method /= PHONON_GAMMA) then
!           write(nfmode,'(1x,"iq=",i5," q=(",f15.10,",",f15.10,",",f15.10,") (",f15.10,",",f15.10,",",f15.10,")")') &
                 & iq,qvin(iq,1:3),qvec(iq,1:3)
         end if
         do i=1,nmodes
***************
*** 2288,2294 ****
        if(mype /= 0) return
        call m_Files_open_nfphdos()
  
!       write(nfphdos,'("# Index Omega(mHa) Omega(eV) Omega(cm-1) DOS(States/Ha) DOS(States/eV) DOS(States/cm-1) IntDOS(States)")') 
        do ie = 0, newin
           eev = e(ie)*ha2ev/1.0d3  ! Hartree
           ecminv = eev*ev2cminv
--- 2288,2294 ----
        if(mype /= 0) return
        call m_Files_open_nfphdos()
  
!       write(nfphdos,'("# Index Omega(mHa) Omega(eV) Omega(cm-1) DOS(States/mHa) DOS(States/eV) DOS(States/cm-1) IntDOS(States)")') 
        do ie = 0, newin
           eev = e(ie)*ha2ev/1.0d3  ! Hartree
           ecminv = eev*ev2cminv
diff -rcN phase0_2019.01/src_phase_3d/m_PseudoPotential.F90 phase0_2019.02/src_phase_3d/m_PseudoPotential.F90
*** phase0_2019.01/src_phase_3d/m_PseudoPotential.F90	2019-03-19 15:15:39.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_PseudoPotential.F90	2019-07-26 09:09:55.094188535 +0900
***************
*** 9438,9443 ****
--- 9438,9444 ----
  !         do iopr=1,nopr
            do iopr=1,nopr+af !ASMS
               ja=napt(ia,iopr)
+              if(ja>natm) ja = ja-natm
               nrorb(ilmta,iopr) = nylm(isph,iopr)
               do mm=1,nylm(isph,iopr)
  ! debug
diff -rcN phase0_2019.01/src_phase_3d/m_Stress.F90 phase0_2019.02/src_phase_3d/m_Stress.F90
*** phase0_2019.01/src_phase_3d/m_Stress.F90	2019-03-19 15:15:39.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_Stress.F90	2019-07-26 09:09:55.084188167 +0900
***************
*** 362,372 ****
  !!$ASASASASAS
  !!$        call G_dot_R(natm,ia,pos,kgp,nbmx,ngabc,zfcos,zfsin)
          if ( kv3/nspin == 1 ) then
-           call G_dot_R(ia,mapmode,1)
            mapmode = MAPPED
          else
-           call G_dot_R(ia,mapmode)
            mapmode = NOTMAPPED
          endif
  !!$ASASASASAS
          do ispin = 1, nspin, af+1
--- 362,372 ----
  !!$ASASASASAS
  !!$        call G_dot_R(natm,ia,pos,kgp,nbmx,ngabc,zfcos,zfsin)
          if ( kv3/nspin == 1 ) then
            mapmode = MAPPED
+           call G_dot_R(ia,mapmode,1)
          else
            mapmode = NOTMAPPED
+           call G_dot_R(ia,mapmode)
          endif
  !!$ASASASASAS
          do ispin = 1, nspin, af+1
diff -rcN phase0_2019.01/src_phase_3d/m_Total_Energy.F90 phase0_2019.02/src_phase_3d/m_Total_Energy.F90
*** phase0_2019.01/src_phase_3d/m_Total_Energy.F90	2019-03-19 15:15:38.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_Total_Energy.F90	2019-07-26 09:09:55.075187763 +0900
***************
*** 709,724 ****
  !!$                  & , trim(tag_mixing), trim(cdmixing_names_applied(1))
            end if
  81        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.12,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 82        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.11,' EDEL = ',d14.6, ' : SOLVER = ',A, ' : ',A,' = ',A)
! 83        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.10,' EDEL = ',d14.6, ' : SOLVER = ',A, ' : ',A,' = ',A)
! 84        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.9,' EDEL = ',d14.6, ' : SOLVER = ',A, ' : ',A,' = ',A)
! 85        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.8,' EDEL = ',d14.6, ' : SOLVER = ',A, ' : ',A,' = ',A)
! 86        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.7,' EDEL = ',d14.6, ' : SOLVER = ',A, ' : ',A,' = ',A)
! 87        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.6,' EDEL = ',d14.6, ' : SOLVER = ',A, ' : ',A,' = ',A)
! 88        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.5,' EDEL = ',d14.6, ' : SOLVER = ',A, ' : ',A,' = ',A)
! 89        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',d20.12,' EDEL = ',d14.6, ' : SOLVER = ',A, ' : ',A,' = ',A)
         else
            write(nfout,600) iteration,etotal,edel
         end if
  
         if(sw_output_xc_seperately==OFF)then
--- 709,731 ----
  !!$                  & , trim(tag_mixing), trim(cdmixing_names_applied(1))
            end if
  81        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.12,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 82        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.11,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 83        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.10,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 84        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.9,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 85        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.8,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 86        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.7,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 87        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.6,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 88        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',f20.5,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
! 89        format(' TOTAL ENERGY FOR',i6,' -TH ITER=',d20.12,' EDEL = ',d14.6, ' : SOLVER = ',A,' + ',A,' : ',A,' = ',A)
         else
            write(nfout,600) iteration,etotal,edel
+ !!$          write(nfout,*) " TOTAL ENERGY FOR ",iteration
+ !!$          call flush(nfout)
+ !!$          write(nfout,*) " -TH ITER", etotal
+ !!$          call flush(nfout)
+ !!$          write(nfout,*) "    edel = ", edel
+ !!$          call flush(nfout)
+ !!$          write(nfout,'("  TOTAL ENERGY FOR",i6," -TH ITER=",f20.8,"   edel = ",d14.6)') iteration, etotal, edel
         end if
  
         if(sw_output_xc_seperately==OFF)then
***************
*** 783,794 ****
  #ifdef __TIMER_SUB__
      call timer_end(749)
  #endif
! 600 FORMAT(' ','TOTAL ENERGY FOR',I6,' -TH ITER=',F20.12,2x,' edel = ',D14.6)
  610 FORMAT(' KI=',F20.12,' HA=',F20.12,' XC=',F20.12,' LO=',F20.12,/ &
!     &      ,' NL=',F20.12,' EW=',F20.12,' PC=',F20.12,' EN=',F20.12)
! 615 FORMAT(' KI=',F20.12,' HA=',F20.12,' EX=',F20.12,' CR=',F20.12,' XC=',F20.12,' LO=',F20.12,/ &
!     &      ,' NL=',F20.12,' EW=',F20.12,' PC=',F20.12,' EN=',F20.12)
! 620 FORMAT(' ','PHYSICALLY CORRECT ENERGY = ',F20.12)
  630 FORMAT(' VD=',F15.7,' VE=',F15.7,' ED=',F15.7)
  640 FORMAT(" HE=",F15.7," HP=",F15.7)
  650 FORMAT(" VEXX=",F20.12," EEXX=",F20.12," EXX-VEXX=",F20.12)
--- 790,801 ----
  #ifdef __TIMER_SUB__
      call timer_end(749)
  #endif
! 600 FORMAT(' ','TOTAL ENERGY FOR',I6,' -TH ITER=',F0.12,2x,' edel = ',D14.6)
  610 FORMAT(' KI=',F20.12,' HA=',F20.12,' XC=',F20.12,' LO=',F20.12,/ &
!     &      ,' NL=',F20.12,' EW=',F20.12,' PC=',d20.12,' EN=',F20.12)
! 615 FORMAT(' KI=',F20.12,' HA=',F20.12,' EX=',d20.12,' CR=',F20.12,' XC=',F20.12,' LO=',F20.12,/ &
!     &      ,' NL=',F20.12,' EW=',F20.12,' PC=',d20.12,' EN=',F20.12)
! 620 FORMAT(' ','PHYSICALLY CORRECT ENERGY = ',d20.12)
  630 FORMAT(' VD=',F15.7,' VE=',F15.7,' ED=',F15.7)
  640 FORMAT(" HE=",F15.7," HP=",F15.7)
  650 FORMAT(" VEXX=",F20.12," EEXX=",F20.12," EXX-VEXX=",F20.12)
diff -rcN phase0_2019.01/src_phase_3d/m_vdWDF.F90 phase0_2019.02/src_phase_3d/m_vdWDF.F90
*** phase0_2019.01/src_phase_3d/m_vdWDF.F90	2019-03-19 15:15:38.000000000 +0900
--- phase0_2019.02/src_phase_3d/m_vdWDF.F90	2019-07-26 09:09:55.074187702 +0900
***************
*** 480,485 ****
--- 480,486 ----
  
    subroutine alloc_vdw()
      integer :: i,cix,ciy,ciz,itmp,c1,c2,c3,ca,cb,cc
+     integer :: imin, imax, maxcount
      real(kind=DP) :: Ta,Tb,Tc, vec(3), bb(3,3)
  
      allocate(imp_fftcd(nabc));imp_fftcd=0
***************
*** 505,513 ****
         allocate(mp_fftcd(nfft_div_size));mp_fftcd=0
         mp_fftcd(:) = mp_fftcd_z(:)
      end if
!     do i=1,nfft_div_size
!        imp_fftcd(mp_fftcd(i)) = i
!     enddo
  #endif
      allocate(ica(nfft_div_size));ica=0
      allocate(icb(nfft_div_size));icb=0
--- 506,537 ----
         allocate(mp_fftcd(nfft_div_size));mp_fftcd=0
         mp_fftcd(:) = mp_fftcd_z(:)
      end if
! !!$ T.Yamasaki 2019.06.12
! !!$    write(nfout,'(" ** alloc_vdw ** nabc = ",i10)') nabc
! !!$    write(nfout,'(" **   sw_fft_xzy = ",i8)') sw_fft_xzy
! !!$    write(nfout,'(" **   nfft_div_size = ",i5)') nfft_div_size
!     imin = 10000
!     imax = -1
!     maxcount = 0
!     do i = 1, nfft_div_size
!        if(imax < mp_fftcd(i)) imax = mp_fftcd(i)
!        if(imin > mp_fftcd(i)) imin = mp_fftcd(i)
!        if(mp_fftcd(i) > nabc) maxcount = maxcount+1
!     end do
! !!$    write(nfout,'(" **   mp_fftcd = [",i8," : ",i8,"], counter of exceeded = ",i8)') imin, imax,maxcount
! !!$    call flush(nfout)
!     
!     if(maxcount > 0) then
!        do i=1,nfft_div_size
!           if(mp_fftcd(i) > nabc) cycle
!           imp_fftcd(mp_fftcd(i)) = i
!        enddo
!     else
!        do i=1,nfft_div_size
!           imp_fftcd(mp_fftcd(i)) = i
!        enddo
!     end if
! !!$ <==
  #endif
      allocate(ica(nfft_div_size));ica=0
      allocate(icb(nfft_div_size));icb=0
***************
*** 685,690 ****
--- 709,715 ----
    subroutine build_theta()
      use progress_bar, only : reset_progress,progress,set_end
      integer :: cqa,ierr
+ !!$    integer :: imax, imin
      integer :: id_sname = -1
      real(kind=DP),allocatable,dimension(:) :: tmpdr,tmpddr
  
***************
*** 707,715 ****
--- 732,757 ----
         allocate(tmpdr(1))
         allocate(tmpddr(1))
      endif
+ !!$    write(nfout,'(" ** build_theta ** ista_nq,iend_nq = ",2i10)') ista_nq,iend_nq
+ !!$    imin = 1000
+ !!$    imax = -1
+ !!$    do cqa = ista_nq,iend_nq
+ !!$       if(imax<map_z_nq(cqa)) imax = map_z_nq(cqa)
+ !!$       if(imin>map_z_nq(cqa)) imin = map_z_nq(cqa)
+ !!$    end do
+ !!$    write(nfout,'(" **   map_z_nq = [",i8,", : ",i8,"]")') imin, imax
+ !!$    write(nfout,'(" **   nfft_div_size = ",i5)') nfft_div_size
+ !!$    if(sw_save_memory_vdw) then
+ !!$       write(nfout,'(" ** np_nq = ",i8)') np_nq
+ !!$    else
+ !!$       write(nfout,'(" ** nq0 = ",i8)') nq0
+ !!$    end if
+ 
      Do cqa = ista_nq,iend_nq
         if(lprog) call progress()
         Call theta_ab(nfft_div_size,cqa,nq0,q0min,q0max,rho,grad,rhomin,theta_R,tmpdr,tmpddr)
+ !!$       write(nfout,'(" ** build_theta na,nb,nc = ",3i8)') na,nb,nc
+ !!$       write(nfout,'("               size of theta_G, theta_R = ",2i8)') nfft_div_size,nfft_div_size
         Call RtoG(na,nb,nc,theta_R,theta_G)
         if(sw_save_memory_vdw)then
            theta_G_ab(map_z_nq(cqa),:) = theta_G(:)
***************
*** 1507,1518 ****
--- 1549,1564 ----
      allocate(temp_R(nabc));temp_R=(0.d0,0.d0)
      allocate(temp_G(nabc));temp_G=(0.d0,0.d0)
      !***** FFT **************************************************
+ !!$    goto 1001
      do ca=1,nfft_div_size
+        if(mp_fftcd(ca)>nabc) cycle
         temp_R(mp_fftcd(ca)) = DCMPLX(theta_R(ca))
      enddo
      call mpi_allreduce(MPI_IN_PLACE,temp_R,nabc,mpi_double_complex,mpi_sum,mpi_ke_world,ierr)
  
      allocate(afft(nfftp_nonpara));afft=0.d0
+ 
+ 
      call map_rho_to_afft(temp_R,afft)
      call m_FFT_CD0(nfout,afft,DIRECT)
      call map_afft_to_rho(afft,temp_G)
***************
*** 1548,1555 ****
  !    enddo
      !***** END of FFT ******************************************
  
!     deallocate(temp_R);   deallocate(temp_G)
!     deallocate(afft)
    End subroutine RtoG
    !** End SUBROUTINE RtoG ***********************************************************************
  
--- 1594,1603 ----
  !    enddo
      !***** END of FFT ******************************************
  
! 1001 continue
!     if(allocated(temp_R)) deallocate(temp_R)
!     if(allocated(temp_G)) deallocate(temp_G)
!     if(allocated(afft))   deallocate(afft)
    End subroutine RtoG
    !** End SUBROUTINE RtoG ***********************************************************************