!{\src2tex{textfont=tt}}
!!****f* ABINIT/wffile
!! NAME
!! wffile
!!
!! FUNCTION
!! Part of cut3d that gives the wavefunction for one kpt,one band
!! and one spin polarisation in real space.  The output depends on
!! the chosen option.
!!
!! COPYRIGHT
!! Copyright (C) 2001-2005 ABINIT group (JFB, MCote, MVer)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~ABINIT/Infos/copyright
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~ABINIT/Infos/contributors .
!!
!! INPUTS
!! Needs an unformatted wave function from abinit.
!! exchn2n3=if 1, n2 and n3 are exchanged
!! headform=format of the wf file
!!
!! natom = number of atoms in the unit cell
!! nr1,nr2,nr3 = grid size (nr1 x nr2 x nr3 = filrho dimension)
!! ntypat = number of atom type
!! ucvol = unit cell volume (> 0)
!! densfileformat = flag for the format of the density file:
!!       0 = ASCII
!!       1 = binary
!! denval = density value exported by interpol, to be wrote in the output file
!! filrho = name of the density file (ASCII or binary)
!! filtau = name of the atomic position file (Xmol format)
!! tau = cartesian coordinates
!! acell = unit cell parameters
!! rprim = orientation of the unit cell axes
!!
!! OUTPUT
!! Depends on the option chosen.
!! It is the wave function for the k point, band and spin polarisation
!! chosen.  It can be written in different ways. The option are describe
!! with the option list.  It is possible to output a Data Explorer file.
!!
!! PARENTS
!!      cut3d
!!
!! CHILDREN
!!      date_and_time,dens_in_sph,fourwf,getkpgnorm,getph,handle_ncerr,hdr_skip
!!      init_bess_spl,initylmg,int2char,kpgio,leave_new,metric,ph1d3d,recip_ylm
!!      rwwf,sort_dp,sphereboundary,splint,wrtout,xredxcart
!!
!! SOURCE

subroutine wffile(ecut,exchn2n3,headform,istwfk,kpt,natom,nband,nkpt,npwarr,&
     &nr1,nr2,nr3,nspinor,nsppol,ntypat,rprimd,tau,typat,znucl)

 use defs_basis
  use defs_datatypes
  use defs_interfaces
 
 implicit none
 
  !no_abirules
  ! For NetCDF
# if defined NETCDF
  use netcdf
# endif
!Arguments -----------------------------------
!scalars
 integer,intent(in) :: exchn2n3,headform,natom,nkpt,nr1,nr2,nr3,nspinor,nsppol
 integer,intent(in) :: ntypat
 real(dp),intent(in) :: ecut
!arrays
 integer,intent(in) :: istwfk(nkpt),nband(nkpt),npwarr(nkpt),typat(natom)
 real(dp),intent(in) :: kpt(3,nkpt),rprimd(3,3),znucl(ntypat)
 real(dp),intent(inout) :: tau(3,natom)
 
!Local variables-------------------------------
  character(len=*), parameter :: INPUTfile='cut.in'
  character(len=fnlen) :: ylmnam=""
!scalars
 integer,save :: tim_fourwf=0,tim_rwwf=0
 integer :: cband,cgshift,ckpt,cplex,cspinor,csppol,formeig,gridshift1
 integer :: gridshift2,gridshift3,ia,iatom,iband,icg,ichoice,ierr,ifile,ii,ii1
 integer :: ii2,ii3,ikpt,ilang,init_prefact,ioffkg,ios,iout,iprompt,ipw,ir1,ir2
 integer :: ir3,ispden,isppol,ivect,ixfh,ixint,mband,mbess,mcg,mgfft,mkmem
 integer :: mlang,mpw,n4,n5,n6,nband_disk,nfit,npw_k,npwout,nradint,nspden,nxfh
 integer :: oldcband,oldckpt,oldcspinor,oldcsppol,option,prtsphere,select_exit
 integer :: unkg,unout=12,unylm=0
 real(dp) :: arg,bessargmax,bessint_delta,kpgmax,ratsph,tmpi,tmpr,ucvol,weight
 real(dp) :: xnow,ynow,znow
 character(len=1) :: outputchar
 character(len=10) :: string
 character(len=4) :: mode_paral
 character(len=500) :: message
 character(len=fnlen) :: kgnam,output,output1
 type(MPI_type) :: mpi_enreg
 type(wffile_type) :: wff
!arrays
 integer :: atindx(natom),iatsph(natom),ngfft(18)
 integer,allocatable :: gbound(:,:),iindex(:),kg(:,:),kg_dum(:,:),kg_k(:,:)
 integer,allocatable :: npwarr1(:),npwtot1(:)
 real(dp) :: acell(3),cmax(natom),gmet(3,3),gprimd(3,3),oldkpt(3)
 real(dp) :: phkxred(2,natom),rmet(3,3),shift_tau(3),tau2(3,natom)
 real(dp) :: xred(3,natom),ylmgr_dum(1)
 real(dp),allocatable :: bess_fit(:,:,:),bess_spl(:,:),bess_spl_der(:,:)
 real(dp),allocatable :: cg(:,:),cgcband(:,:),denpot(:,:,:),eigen(:)
 real(dp),allocatable :: fofgin(:,:),fofgout(:,:),fofr(:,:,:,:),kpgnorm(:)
 real(dp),allocatable :: occ1(:),ph1d(:,:),ph3d(:,:,:),rint(:),spl_bessint(:,:)
 real(dp),allocatable :: sum_1atom_1ll(:,:),x_bess(:),xfhist(:,:,:,:),xfit(:)
 real(dp),allocatable :: yfit(:),ylm(:,:),ylm_k(:,:)
 character(len=fnlen),allocatable :: filename(:)
  !no_abirules
  !For NetCDF********************************************************
# if defined NETCDF
  integer :: ncid, status, gridsize1DimID, gridsize2DimID, gridsize3DimID
  integer :: latDimID, nbatomDimID, imagwavefunVarID,realwavefunVarID, latticevecVarID, originVarID,grid1VarID,grid2VarID,grid3VarID
  integer :: atomposiVarID, atomicnumVarID, titlechoice, posDimID,kpointVarID 
  character*500 :: filetitle 
  integer :: dd,mm,yyyy, igrid
  integer :: values(8)
  character(len=5) :: strzone
  character(len=8) :: strdat
  character(len=10) :: strtime
  character(len=3), parameter :: monnam(12)=(/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'/)
  character(len=11) ::stridate
  real :: originatt(3,3), gridwavefun1(3,2), gridwavefun2(3,2), gridwavefun3(3,2)
  real,allocatable :: partwf(:,:,:)
  real :: kptvar(3)
# endif
 
! ***********************************************************************

! BEGIN EXECUTABLE SECTION

  mpi_enreg%paralbd=0
  mpi_enreg%nproc_fft=1
  mpi_enreg%me=0
  mpi_enreg%me_fft=0
  mpi_enreg%paral_fft=0
  mpi_enreg%paral_compil_kpt=0
  mpi_enreg%paral_compil_fft=0
  mpi_enreg%fft_option_lob=1

  formeig=0
  oldckpt=0
  oldcband=0
  oldcsppol=0
  oldcspinor=0

  iout=-1
  call metric(gmet,gprimd,iout,rmet,rprimd,ucvol)

  ! get xred
  call xredxcart(natom,-1,rprimd,tau,xred)

  do iatom=1,natom
     iatsph(iatom) = iatom
     atindx(iatom) = iatom
  end do

  ! max ang mom + 1
  mlang = 5

  allocate (kg_dum(3,0))

  allocate (ph1d(2,(2*nr1+1+2*nr2+1+2*nr3+1)*natom))
  call getph(atindx,natom,nr1,nr2,nr3,ph1d,xred)

  do
     !Get k-point, band and spin polarisation for the output
     if(nkpt/=1)then
        write(*,*)
        write(*,'(a,i4,a)') ' For which k-points? (1 to ',nkpt,')'
        read(5,*)ckpt
        !Check if kpt exist
        if(ckpt<1 .or. ckpt>nkpt) then
           write(6,*) 'Invalid k-point ',ckpt
           stop
        end if
     else
        ckpt=nkpt
     end if
     write(*,*) ' => Your k-point is : ',ckpt
     write(*,*) '    Reduced coordinates: ', &
        & kpt(1,ckpt), '  ', kpt(2, ckpt), '  ',        &
        & kpt(3, ckpt)
     write(*,*)

     if(nband(ckpt)/=1)then
        write(*,*)
        write(*,'(a,i5,a)') ' For which band ? (1 to ',nband(ckpt),')'
        read(5,*)cband
        !Check if band number exist

        if(cband<1 .or. cband>nband(ckpt)) then
           write(6,*) 'Invalid band number',cband
           stop
        end if
     else
        cband=nband(ckpt)
     end if
     write(*,*) ' => Your band number is : ',(cband)
     write(*,*)

     if(nsppol/=1)then
        write(*,*)
        write(*,*) ' For which spin polarisation ?'
        read(5,*)csppol
        !Check if spin polarisation exist
        if(csppol<1 .or. csppol>nsppol) then
           write(6,*)'Invalid spin polarisation ',csppol
           stop
        end if
     else
        csppol=1
     end if

     write(*,*) ' => Your spin polarisation number is : ',(csppol)
     write(*,*)

     if(nspinor/=1) then
        write (*,*) ' nspinor = ', nspinor
        write(*,*)
        write(*,*) ' For which spinor component ?'
        read(5,*) cspinor
        !Check if spin polarisation exist
        if(cspinor<1 .or. cspinor>nspinor) then
           write(6,*)'Invalid spinor index ',cspinor
           stop
        end if
        write(*,*) ' => Your spinor component is : ',(cspinor)
        write(*,*)
     else
        cspinor=1
     end if

     !Reading of the data if the value of ckpt and csppol are
     !different from oldckpt and oldcsppol
     !    formeig=0 gstate calculation
     ! else formeig=1 for response function calculation
     if(csppol/=oldcsppol .or. ckpt/=oldckpt)then
        mband=maxval(nband)
        mpw=maxval(npwarr)
        mcg=mpw*nspinor*mband
        if (allocated (cg))  deallocate (cg,eigen,occ1)
        allocate(cg(2,mcg),eigen((2*mband)**formeig*mband),&
             & occ1(mband))
        !Rewind the file and skip the header
        wff%unwff=19
        wff%accesswff=0
        call hdr_skip(wff,ierr)
        !Should use xdefineOff if MPI I/O
        !Iteration on nsppol and kpt to skip the unwanted datas (option -1)
        do isppol=1,csppol
           do ikpt=1,nkpt
              if(isppol==csppol .and. ikpt==ckpt)then
                 option=1
              else
                 option=-1
              end if
              call rwwf(cg,eigen,formeig,headform,0,ikpt,isppol,kg_dum,&
                   &          mband,mcg,nband(ikpt),nband_disk,&
                   &          npwarr(ikpt),nspinor,occ1,option,0,tim_rwwf,wff)

              !         In case one read the last wf, the history will be output too
              if(ikpt==nkpt .and. isppol==nsppol)then
                 if(cband==nband(nkpt))then
                    read(unit=19,iostat=ios)nxfh
                    if(ios>0)then
                       write(message, '(a,a,a,a,a,a)' )ch10,&
                            &             ' wffile : BUG -',ch10,&
                            &             '  An error occurred reading the input wavefunction file,',ch10,&
                            &             '  history record.'
                       call wrtout(6,message,'COLL')
                       call leave_new('COLL')
                    else if(ios==0)then
                       write(message, '(a,a,i4,a)' )ch10,&
                            &             ' wffile : reading',nxfh,&
                            &             ' (x,f) history pairs from input wf file.'
                       call wrtout(6,message,'COLL')
                    end if
                    if(nxfh>=1)then
                       allocate(xfhist(3,natom+4,2,nxfh))
                       do ixfh=1,nxfh
                          read(19)xfhist(:,:,:,ixfh)
                          write(message, '(a,a,i6,a)')ch10,&
                               &              ' History step number ',ixfh,&
                               &              ' , atom number, xcart(1:3), fcart(1:3) ='
                          call wrtout(6,message,'COLL')
                          do iatom=1,natom
                             write(6,'(i5,6es16.6)')iatom,xfhist(:,iatom,:,ixfh)
                          end do
                       end do
                       deallocate(xfhist)
                    end if ! nxfh>1
                 end if ! cband==nband(nkpt)
              end if ! ikpt==nkpt .and. isppol==nsppol

              if(option==1)exit  ! When the target wf has been read, exit the wf file reading
           end do
           if(option==1)exit
        end do

     end if

     if (csppol/=oldcsppol .or. ckpt/=oldckpt .or. &
          & cband/=oldcband .or. cspinor/=oldcspinor ) then
        !The data of ckpt,cnsspol are in cg
        !Now we have to do the Fourier Transform of the datas

        ngfft(1)=nr1
        ngfft(2)=nr2
        ngfft(3)=nr3
        !ngfft(4) and ngfft(5) can not be even (see getng.f)
        if (mod(nr1,2)==0)then
           ngfft(4)=nr1+1
        else
           ngfft(4)=nr1
        end if
        if (mod(nr2,2)==0)then
           ngfft(5)=nr2+1
        else
           ngfft(5)=nr2
        end if
        ngfft(6)=nr3
        !       XG 020829 : 112 does not work yet for all istwfk values
        ngfft(7)=111
        ngfft(8)=256
        ngfft(9)=0
        ngfft(10)=1
        ngfft(11)=0
        ngfft(12)=ngfft(2)
        ngfft(13)=ngfft(3)
        ngfft(14)=0

        !if iout<0, the output of metric will not be print
        mode_paral='PERS'
        mkmem=nkpt
        mgfft=maxval(ngfft(1:3))
        allocate(npwarr1(nkpt),kg(3,mpw*mkmem),npwtot1(nkpt))

        !       Create positions index for pw
        call kpgio(ecut,exchn2n3,gmet,istwfk,kg,kgnam,kpt,mkmem,nband,nkpt,&
             &         mode_paral,mpi_enreg,mpw,npwarr1,npwtot1,nsppol,unkg)

        allocate(ylm(mpw*nkpt,mlang*mlang))
        call initylmg(gprimd,kg,kpt,mkmem,mpi_enreg,&
             &            mlang,mpw,nband,nkpt,npwarr,nsppol,0,rprimd,unkg,unylm,ylm,ylmgr_dum)

        ioffkg=0
        do ikpt=1,ckpt-1
           ioffkg=ioffkg+npwarr1(ikpt)
        end do
        npw_k=npwarr(ckpt)
        allocate(gbound(2*mgfft+8,2),kg_k(3,npw_k))
        kg_k(:,1:npw_k)=kg(:,1+ioffkg:npw_k+ioffkg)

        !       Compute the norms of the k+G vectors
        allocate(kpgnorm(npw_k))
        call getkpgnorm(gprimd,kpt(:,ckpt),kg_k,kpgnorm,npw_k)

        call sphereboundary(gbound,istwfk(ckpt),kg_k,mgfft,npw_k)
        !Do the Fourier Transform
        n4=ngfft(4)
        n5=ngfft(5)
        n6=ngfft(6)
        !        cplex=0
        cplex=1
        !Complex can be set to 0 with this option(0) of fourwf
        !The shift is to get the good band values
        cgshift=(cband-1)*npw_k*nspinor + (cspinor-1)*npw_k
        allocate(cgcband(2,npw_k))
        cgcband(:,1:npw_k)=cg(:,cgshift+1:cgshift+npw_k)

        !       Fix the phase of cgcband, for portability reasons
        !       call fxphas(cgcband,cgcband,0,npw_k,1,npw_k,0)

        allocate(denpot(cplex*n4,n5,n6))
        allocate(fofgout(2,npw_k))
        allocate(fofr(2,n4,n5,n6))
        call fourwf(cplex,denpot,cgcband,fofgout,fofr,&
             &  gbound,gbound,&
             &  istwfk(ckpt),kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,npw_k,&
             &  npw_k,n4,n5,n6,0,tim_fourwf,weight)

        !       Analyse wavefunction inside atomic sphere

        write (6, '(a)' ) ' Do you want the atomic analysis for this state : '
        write (6, '(a,2i5,a,a)' ) ' (kpt,band)= (',ckpt,cband,')? '
        write (6, '(a)' ) ' If yes, enter the radius of the atomic spheres, in bohr '
        write (6, '(a)' ) ' If no, enter 0 '
        read (*,*) ratsph
        write(6, '(a,f16.8,a)' ) ' You entered ratsph=',ratsph,' Bohr '

        if (ratsph >= tol10) then

           write(6, '(3a)' ) ch10,' Atomic sphere analysis ',ch10

           !        Init bessel function integral for recip_ylm
           !        max ang mom + 1
           mlang = 5
           bessint_delta = 0.1_dp
           kpgmax = sqrt(ecut)
           bessargmax = ratsph*two_pi*kpgmax
           mbess = int (bessargmax / bessint_delta) + 1
           bessargmax = bessint_delta*mbess

           !        Intervals in radial integration
           nradint = mbess

           write (6, '(a,2es16.6,i6)' ) &
                &         ' wffile : kpgmax, bessargmax, nradint = ', kpgmax, bessargmax,nradint

           !        Initialize general Bessel function array on uniform grid
           !        x_bess, from 0 to (2 \pi |k+G|_{max} |r_{max}|)
           allocate (bess_spl(mbess,mlang))
           allocate (bess_spl_der(mbess,mlang))
           allocate (x_bess(nradint),rint(nradint))

           call init_bess_spl(mbess,bessargmax,bessint_delta,mlang,&
                &         bess_spl,bess_spl_der,x_bess)
           !DEBUG
           ! write (*,*) 'wffile : after init_bess_spl :'
           ! write (*,'(6F12.5)') bess_spl
           !ENDDEBUG

           allocate (bess_fit(mpw,nradint,mlang),xfit(npw_k),yfit(npw_k))
           allocate (iindex(npw_k))
           nfit = npw_k

           do ixint=1,nradint
              rint(ixint) = (ixint-1)*ratsph / (nradint-1)

              do ipw=1,npw_k
                 xfit(ipw) = two_pi * kpgnorm(ipw) * rint(ixint)
                 iindex(ipw) = ipw
              end do
              call sort_dp (npw_k,xfit,iindex,tol14)
              do ilang=1,mlang
                 call splint(mbess,x_bess,bess_spl(:,ilang),bess_spl_der(:,ilang),&
                      &                  nfit,xfit,yfit)
                 !          Re-order results for different G vectors
                 do ipw=1,npw_k
                    bess_fit(iindex(ipw),ixint,ilang) = yfit(ipw)
                 end do
              end do ! ipw
           end do ! ixint

           !        Construct phases ph3d for all G vectors in present sphere
           !        make phkred for all atoms

           do ia=1,natom
              iatom=atindx(ia)
              arg=two_pi*( kpt(1,ckpt)*xred(1,ia) &
                   &                    + kpt(2,ckpt)*xred(2,ia) &
                   &                    + kpt(3,ckpt)*xred(3,ia))
              phkxred(1,iatom)=cos(arg)
              phkxred(2,iatom)=sin(arg)
           end do

           allocate (ph3d(2,npw_k,natom))
           !        Get full phases for the following
           !        write (*,*) 'nr1nr2nr3 ',nr1,nr2,nr3
           call ph1d3d(1,natom,kg_k,kpt(:,ckpt),natom,natom,npw_k,nr1,nr2,nr3,&
                &            phkxred,ph1d,ph3d)
           !        phases exp (2 pi i (k+G).x_tau) are now in ph3d

           allocate (ylm_k(mpw,mlang*mlang))
           do ilang=1,mlang*mlang
              do ipw=1,npw_k
                 ylm_k(ipw,ilang) = ylm(ioffkg+ipw,ilang)
              end do
           end do
           allocate(sum_1atom_1ll(mlang,natom))
           prtsphere=1
           call recip_ylm (bessargmax,bess_fit,cgcband,iatsph,&
                &                istwfk(ckpt),kg_k,kpgnorm,&
                &                nradint,mgfft,mlang,mpi_enreg,mpw,natom,natom,ngfft,npw_k,&
                &                ntypat,ph3d,prtsphere,rint,&
                &                ratsph,rprimd,sum_1atom_1ll,typat,ucvol,ylm_k,znucl)

           call dens_in_sph(cmax,cgcband,gmet,istwfk(ckpt),&
                &           kg_k,natom,ngfft,mpi_enreg,npw_k,ph1d,ratsph,ucvol)

           write(6, '(a)' )' Charge in the sphere around each atom '
           do iatom=1,natom
              write(6, '(a,i4,a,f14.8)' ) ' Atom number ',iatom,' :  charge =',cmax(iatom)
           end do

           deallocate(sum_1atom_1ll)
           deallocate(ylm_k,ph3d,iindex,yfit,xfit,bess_fit)
           deallocate (bess_spl)
           deallocate (bess_spl_der)
           deallocate (x_bess,rint)
        end if ! ratsph < 0     = end if for atomic sphere analysis

        deallocate(cgcband,fofgout,denpot,gbound,kg_k)
        deallocate(npwarr1,kg,npwtot1,kpgnorm,ylm)

     end if


     write(6,*)
     write(6,*) ' 3D wave function was read. ',&
          &               'Ready for further treatment.'
     write(6,*)
     write(6,*) '============================',&
          &               '==============================='
     write(6,*)

     !------------------------------------------------------------------------

     ! At this moment all the input is done
     ! The code knows the geometry of the system,
     ! and the data file.


     select_exit = 0
     do while (select_exit == 0)
        write(*,*) ' What is your choice ? Type:'
        write(*,*) '  1 => 3D formatted real and imaginary data'
        write(*,*) '       (output the bare 3D data - two column,R,I)'
        write(*,*) '  2 => 3D formatted real data'
        write(*,*) '       (output the bare 3D data - one column)'
        write(*,*) '  3 => 3D formatted imaginary data'
        write(*,*) '       (output the bare 3D data - one column)'
        write(*,*) '  4 => 3D indexed real and imaginary data'
        write(*,*) '       (3D data, preceeded by 3D index)'
        write(*,*) '  5 => 3D indexed real data'
        write(*,*) '       (bare 3D data, preceeded by 3D index)'
        write(*,*) '  6 => 3D indexed imaginary data'
        write(*,*) '       (bare 3D data, preceeded by 3D index)'
        write(*,*) '  7 => 3D Data Explorer formatted data '
        write(*,*) '       (Real file and Imaginary file)'
        write(*,*) '  8 => 3D Data Explorer formatted data '
        write(*,*) '       (Only the Real file)'
        write(*,*) '  9 => 3D Data Explorer formatted data '
        write(*,*) '       (Only the Imaginary file)'
        write(*,*) ' 10 => 3D Data Explorer formatted data and position files'
        write(*,*) ' 11 => XCrysden formatted data and position files'
        write(*,*) ' 12 => NetCDF data and position file'
        write(*,*) ' 13 => XCrysden/VENUS wavefunction data'
        write(*,*) ' 14 => exit to k-point / band / spin-pol loop'
        read(*,*) ichoice
        write(*, '(a,a,i2,a)' ) ch10,' Your choice is ',ichoice,char(10)

        if (ichoice>0 .and. ichoice<14)then
           write(*,*) ch10,'  Enter the root of an output file:'
           read(*,*) output1
           write(*,*) '  The root of your file is : ',trim(output1)
           output=trim(output1)//'_k'
           call int2char(ckpt,string)
           output=trim(output)//trim(string)//'_b'
           call int2char(cband,string)
           output=trim(output)//trim(string)//'_s'
           call int2char(csppol,string)
           output=trim(output)//trim(string)
           write(*,*) '  The corresponding filename is : ',trim(output)
        end if

        select case(ichoice)

        case(1)            ! data R,I
           write(*,*)
           write(*,*) 'Give 1 file of 3D formatted real and imaginary data'
           write(*,*) 'The first column is the real data'
           write(*,*) 'The second column is the imaginary data'
           write(*,*)
           open(unit=unout,file=output,status='replace',form='formatted')
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    write(unout,'(2f20.16)') fofr(:,ir1,ir2,ir3)
                 end do
              end do
           end do
           close(unout)
           exit

        case(2)            ! data R
           write(*,*)
           write(*,*) 'Give 1 file of 3D formatted real data'
           write(*,*) 'The only column is the real data'
           write(*,*)
           open(unit=unout,file=output,status='replace',form='formatted')
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    write(unout,'(f20.16)') fofr(1,ir1,ir2,ir3)
                 end do
              end do
           end do
           close(unout)
           exit

        case(3)            ! data I
           write(*,*)
           write(*,*) 'Give 1 file of 3D formatted real data'
           write(*,*) 'The only column is the imaginary data'
           write(*,*)
           open(unit=unout,file=output,status='replace',form='formatted')
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    write(unout,'(f20.16)') fofr(2,ir1,ir2,ir3)
                 end do
              end do
           end do
           close(unout)
           exit

        case(4)            ! coord(x,y,z) data R,I
           write(*,*)
           write(*,*) 'Give 1 file of 3D formatted data'
           write(*,*) 'The first three columns are the x,y,z positions(Angstrom)'
           write(*,*) 'The fourth column is the real data'
           write(*,*) 'The fifth column is the imaginary data'
           write(*,*)
           open(unit=unout,file=output,status='replace',form='formatted')
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    xnow = rprimd(1,1)*(ir1-1)/nr1 + rprimd(1,2)*(ir2-1)/nr2 + rprimd(1,3)*(ir3-1)/nr3
                    ynow = rprimd(2,1)*(ir1-1)/nr1 + rprimd(2,2)*(ir2-1)/nr2 + rprimd(2,3)*(ir3-1)/nr3
                    znow = rprimd(3,1)*(ir1-1)/nr1 + rprimd(3,2)*(ir2-1)/nr2 + rprimd(3,3)*(ir3-1)/nr3
                    write(unout,'(3f16.10,2f20.16)') Bohr_Ang*xnow, Bohr_Ang*ynow, Bohr_Ang*znow,fofr(:,ir1,ir2,ir3)
                 end do
              end do
           end do
           close(unout)
           exit

        case(5)            ! coord(x,y,z) data R
           write(*,*)
           write(*,*) 'Give 1 file of 3D formatted data'
           write(*,*) 'The first three columns are the x,y,z positions(Angstrom)'
           write(*,*) 'The fourth column is the real data'
           write(*,*)
           open(unit=unout,file=output,status='replace',form='formatted')
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    xnow = rprimd(1,1)*(ir1-1)/nr1 + rprimd(1,2)*(ir2-1)/nr2 + rprimd(1,3)*(ir3-1)/nr3
                    ynow = rprimd(2,1)*(ir1-1)/nr1 + rprimd(2,2)*(ir2-1)/nr2 + rprimd(2,3)*(ir3-1)/nr3
                    znow = rprimd(3,1)*(ir1-1)/nr1 + rprimd(3,2)*(ir2-1)/nr2 + rprimd(3,3)*(ir3-1)/nr3
                    write(unout,'(3f16.10,2f20.16)') Bohr_Ang*xnow, Bohr_Ang*ynow, Bohr_Ang*znow,fofr(1,ir1,ir2,ir3)
                 end do
              end do
           end do
           close(unout)

           exit
        case(6)            ! coord(x,y,z) data I
           write(*,*)
           write(*,*) 'Give 1 file of 3D formatted data'
           write(*,*) 'The first three columns are the x,y,z positions(Angstrom)'
           write(*,*) 'The fourth column is the imaginary data'
           write(*,*)
           open(unit=unout,file=output,status='replace',form='formatted')
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    xnow = rprimd(1,1)*(ir1-1)/nr1 + rprimd(1,2)*(ir2-1)/nr2 + rprimd(1,3)*(ir3-1)/nr3
                    ynow = rprimd(2,1)*(ir1-1)/nr1 + rprimd(2,2)*(ir2-1)/nr2 + rprimd(2,3)*(ir3-1)/nr3
                    znow = rprimd(3,1)*(ir1-1)/nr1 + rprimd(3,2)*(ir2-1)/nr2 + rprimd(3,3)*(ir3-1)/nr3
                    write(unout,'(3f16.10,2f20.16)') Bohr_Ang*xnow, Bohr_Ang*ynow, Bohr_Ang*znow,fofr(2,ir1,ir2,ir3)
                 end do
              end do
           end do
           close(unout)
           exit

        case(7)            !OpenDX format, data R and data I
           write(*,*)
           write(*,*) 'Give 2 files of 3D formatted data'
           write(*,*) 'The file is ready to be use with OpenDX'
           write(*,*) 'The eigenvalues and occupations numbers are in comments'
           write(*,*)
           allocate(filename(2))
           filename(1)=trim(output)//'Real.dx'
           filename(2)=trim(output)//'Imag.dx'
           write(*,*) '  The name of your files is : '
           write(*,*) trim(filename(1)),'  for the real part,'
           write(*,*) trim(filename(2)),'  for the imaginary part.'
           write(*,*)

           do ifile=1,2
              open(unit=unout,file=filename(ifile),status='replace',form='formatted')
              rewind(unout)
              write(unout,*)'# band,  eigenvalues   and   occupations'
              do iband=1,nband(ckpt)
                 write(unout,'(a,i4,2f20.16)')'#',iband,eigen(iband),occ1(iband)
              end do
              write(unout,'(a,i10,a)')'object "donnees" class array type float rank 0 items',&
                   &nr1*nr2*nr3,' data follows'
              do ir3=1,nr3
                 do ir2=1,nr2
                    do ir1=1,nr1
                       write(unout,'(f20.16)')fofr(ifile,ir1,ir2,ir3)
                    end do
                 end do
              end do

              write(unout,'(a)')'# this is the object defining the grid connections'
              write(unout,'(a,3i5)')'object "gridconnections" class gridconnections counts',&
                   &nr3,nr2,nr1
              write(unout,*)
              write(unout,*)
              write(unout,'(a)')'# this is the object defining the grid'
              write(unout,'(a,3i5)')'object "positions" class gridpositions counts',&
                   &nr3,nr2,nr1

              write(unout,'(a)') 'origin 0 0 0'
              write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,3)/nr3, ii1=1,3)
              write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,2)/nr2, ii1=1,3)
              write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,1)/nr1, ii1=1,3)

              write(unout,'(a)')'# this is the collective object, one for each grid '
              write(unout,'(a)')'object "densite" class field '
              write(unout,'(a)')'component "positions"   value "positions"'
              write(unout,'(a)')'component "connections" value "gridconnections" '
              write(unout,'(a)')'component "data"        value "donnees"'

              close(unit=unout)
           end do
           deallocate(filename)
           exit

        case(8)            !OpenDX format, data R and data I
           write(*,*)
           write(*,*) 'Give 2 files of 3D formatted data'
           write(*,*) 'The file is ready to be use with OpenDX'
           write(*,*) 'The eigenvalues and occupations numbers are in comments'
           write(*,*)
           allocate(filename(1))
           filename(1)=trim(output)//'Real.dx'
           write(*,*) '  The name of your file is : '
           write(*,*) trim(filename(1)),'  for the real part,'
           write(*,*)


           open(unit=unout,file=filename(1),status='replace',form='formatted')
           rewind(unout)
           write(unout,*)'# band,  eigenvalues   and   occupations'
           do iband=1,nband(ckpt)
              write(unout,'(a,i4,2f20.16)')'#',iband,eigen(iband),occ1(iband)
           end do
           write(unout,'(a,i10,a)')'object "donnees" class array type float rank 0 items',&
                &nr1*nr2*nr3,' data follows'
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    write(unout,'(f20.16)')fofr(1,ir1,ir2,ir3)
                 end do
              end do
           end do

           write(unout,'(a)')'# this is the object defining the grid connections'
           write(unout,'(a,3i5)')'object "gridconnections" class gridconnections counts',&
                &nr3,nr2,nr1
           write(unout,*)
           write(unout,*)
           write(unout,'(a)')'# this is the object defining the grid'
           write(unout,'(a,3i5)')'object "positions" class gridpositions counts',&
                &nr3,nr2,nr1

           write(unout,'(a)') 'origin 0 0 0'
           write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,3)/nr3, ii1=1,3)
           write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,2)/nr2, ii1=1,3)
           write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,1)/nr1, ii1=1,3)

           write(unout,'(a)')'# this is the collective object, one for each grid '
           write(unout,'(a)')'object "densite" class field '
           write(unout,'(a)')'component "positions"   value "positions"'
           write(unout,'(a)')'component "connections" value "gridconnections" '
           write(unout,'(a)')'component "data"        value "donnees"'

           close(unit=unout)
           deallocate(filename)
           exit

        case(9)            !OpenDX format, data R and data I
           write(*,*)
           write(*,*) 'Give 2 files of 3D formatted data'
           write(*,*) 'The file is ready to be use with OpenDX'
           write(*,*) 'The eigenvalues and occupations numbers are in comments'
           write(*,*)
           allocate(filename(1))
           filename(1)=trim(output)//'Imag.dx'
           write(*,*) '  The name of your file is : '
           write(*,*) trim(filename(1)),'  for the imaginary part.'
           write(*,*)


           open(unit=unout,file=filename(1),status='replace',form='formatted')
           rewind(unout)
           write(unout,*)'# band,  eigenvalues   and   occupations'
           do iband=1,nband(ckpt)
              write(unout,'(a,i4,2f20.16)')'#',iband,eigen(iband),occ1(iband)
           end do
           write(unout,'(a,i10,a)')'object "donnees" class array type float rank 0 items',&
                &nr1*nr2*nr3,' data follows'
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    write(unout,'(f20.16)')fofr(2,ir1,ir2,ir3)
                 end do
              end do
           end do

           write(unout,'(a)')'# this is the object defining the grid connections'
           write(unout,'(a,3i5)')'object "gridconnections" class gridconnections counts',&
                &nr3,nr2,nr1
           write(unout,*)
           write(unout,*)
           write(unout,'(a)')'# this is the object defining the grid'
           write(unout,'(a,3i5)')'object "positions" class gridpositions counts',&
                &nr3,nr2,nr1

           write(unout,'(a)') 'origin 0 0 0'
           write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,3)/nr3, ii1=1,3)
           write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,2)/nr2, ii1=1,3)
           write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,1)/nr1, ii1=1,3)

           write(unout,'(a)')'# this is the collective object, one for each grid '
           write(unout,'(a)')'object "densite" class field '
           write(unout,'(a)')'component "positions"   value "positions"'
           write(unout,'(a)')'component "connections" value "gridconnections" '
           write(unout,'(a)')'component "data"        value "donnees"'

           close(unit=unout)
           deallocate(filename)
           exit

        case(10)           !OpenDX format, data R and data I, atoms positions, lattice and cell
           write(*,*)
           write(*,*) 'Give 5 files of formatted data'
           write(*,*) 'The files are ready to be use with Data Explorer'
           write(*,*) 'The eigenvalues and occupations numbers are in comments'
           write(*,*) 'of the two data files'
           write(*,*)
           allocate(filename(2))
           filename(1)=trim(output)//'Real.dx'
           filename(2)=trim(output)//'Imag.dx'
           write(*,*) '  The name of your data files is : '
           write(*,*) trim(filename(1)),'  for the real part,'
           write(*,*) trim(filename(2)),'  for the imaginary part.'
           write(*,*)

           do ifile=1,2
              open(unit=unout,file=filename(ifile),status='replace',form='formatted')
              rewind(unout)
              do iband=1,nband(ckpt)
                 write(unout,'(a,2f20.16)')'#', eigen(iband),occ1(iband)
              end do
              write(unout,'(a,i10,a)')'object "donnees" class array type float rank 0 items',&
                   &nr1*nr2*nr3,' data follows'
              do ir3=1,nr3
                 do ir2=1,nr2
                    do ir1=1,nr1
                       write(unout,'(f20.16)')fofr(ifile,ir1,ir2,ir3)
                    end do
                 end do
              end do

              write(unout,'(a)')'# this is the object defining the grid connections'
              write(unout,'(a,3i5)')'object "gridconnections" class gridconnections counts',&
                   &nr3,nr2,nr1
              write(unout,*)
              write(unout,*)
              write(unout,'(a)')'# this is the object defining the grid'
              write(unout,'(a,3i5)')'object "positions" class gridpositions counts',&
                   &nr3,nr2,nr1

              write(unout,'(a)') 'origin 0 0 0'
              write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,3)/nr3, ii1=1,3)
              write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,2)/nr2, ii1=1,3)
              write(unout,'(a,3f16.10)')'delta ',(Bohr_Ang*rprimd(ii1,1)/nr1, ii1=1,3)

              write(unout,'(a)')'# this is the collective object, one for each grid '
              write(unout,'(a)')'object "densite" class field '
              write(unout,'(a)')'component "positions"   value "positions"'
              write(unout,'(a)')'component "connections" value "gridconnections" '
              write(unout,'(a)')'component "data"        value "donnees"'

              close(unit=unout)
           end do
           deallocate(filename)
           !
           !     write LATTICE_VEC.dx file
           !
           allocate(filename(3))
           filename(1)=trim(output1)//'_LATTICE_VEC.dx'
           filename(2)=trim(output1)//'_ATOM_POS.dx'
           filename(3)=trim(output1)//'_UCELL_FRAME.dx'
           write(*,*)
           write(*,*)'Give the lattice file, ', trim(filename(1))
           open(unit=unout,file=filename(1),status='replace')
           write(unout,'("#",/,"#",/,"#    LATTICE VECTOR INFO:",/,"#",/,"#")')
           write(unout,'(a)') 'object "lattices" class array type float rank 1 shape 3 items 3 data follows'
           do ivect=1,3
              write(unout,'(3f16.10)')  Bohr_Ang*rprimd(1,ivect),Bohr_Ang*rprimd(2,ivect),Bohr_Ang*rprimd(3,ivect)
           end do
           write(unout,'(a,a)') 'object "lattices_location" class array type float ',&
                & 'rank 1 shape 3 items 3 data follows'
           do ivect=1,3
              write(unout,'(3f16.10)')  0d0,0d0,0d0
           end do
           write(unout,'("object   3 class field")')
           write(unout,'(a)') 'component "data" value "lattices"'
           write(unout,'(a)') 'component "positions" value "lattices_location"'
           close(unout)


           !
           !     write ATOM_POS.dx file
           !
           write(*,*)'Give the atoms positions file, ', trim(filename(2))
           open(unit=unout,file=filename(2),status='unknown')
           write(unout,'("#",/,"#",/,"#    BALL AND STICK INFO:",/,"#",/,"#")')
           write(unout,'(a,i5,a)') 'object "atomcoord" array type float rank 1 shape 3 items ',natom,' data follows'
           do iatom=1,natom
              write(unout,'(3f16.10)')  Bohr_Ang*tau(1:3,iatom)
           end do
           !  write(unout,'(a,i5,a)') 'object "data" array type string rank 0 shape 2 items ',&
           !&natom,' data follows'
           write(unout,'(a,i5,a)') 'object "colorcode" array type float rank 0 items ',natom,' data follows'
           do iatom=1,natom
              write(unout,'(f10.4)') znucl(typat(iatom))
           end do
           write(unout,'(a)') 'object "molecule" field'
           write(unout,'(a)') 'component "positions" value "atomcoord"'
           write(unout,'(a)') 'component "data" value "colorcode"'
           close(unout)

           !
           !     write UCELL_FRAME.dx file
           !
           write(*,*)'Give the enveloppe of the cell file, ',trim(filename(3))
           open(unit=unout,file=filename(3),status='unknown')
           write(unout,'("#",/,"#",/,"#    UNIT CELL FRAME INFO:",/,"#",/,"#")')
           write(unout,'(a)')'object 3 class array type int rank 1 shape 2 items 12 data follows'
           write(unout,'(" 0  1",/," 0  2",/," 0  3",/," 1  4",/," 1  5",/," 3  5")')
           write(unout,'(" 3  6",/," 2  6",/," 2  4",/," 7  5",/," 7  6",/," 7  4")')
           write(unout,'(a)') 'attribute "element type" string "lines"'
           write(unout,'("object  4 class array type float rank 1 shape 3 items    8 data follows")')
           write(unout,'("      .00000000      .00000000      .00000000")')
           write(unout,'(3f20.10)') Bohr_Ang*rprimd(:,1)
           write(unout,'(3f20.10)') Bohr_Ang*rprimd(:,2)
           write(unout,'(3f20.10)') Bohr_Ang*rprimd(:,3)
           write(unout,'(3f20.10)') Bohr_Ang*(rprimd(:,1)+rprimd(:,2))
           write(unout,'(3f20.10)') Bohr_Ang*(rprimd(:,1)+rprimd(:,3))
           write(unout,'(3f20.10)') Bohr_Ang*(rprimd(:,2)+rprimd(:,3))
           write(unout,'(3f20.10)') Bohr_Ang*(rprimd(:,1)+rprimd(:,2)+rprimd(:,3))
           write(unout,'("object 5 array type float rank 0 items 12 data follows")')
           do ivect=1,12
              write(unout,'("1.0")')
           end do
           write(unout,'(a)') 'attribute "dep" string "connections"'
           write(unout,'("object 6 class field")')
           write(unout,'(a)') 'component "data" value 5'
           write(unout,'(a)') 'component "positions" value 4'
           write(unout,'(a)') 'component "connections" value 3'
           close(unout)
           deallocate(filename)

           write(*,*)
           exit

        case(11)
           write(*,*)
           write(*,*) 'Give 1 files of formatted data'
           write(*,*) 'The files are ready to be used with XCrysDen'
           write(*,*)
           gridshift1 = 0
           gridshift2 = 0
           gridshift3 = 0
           write(*,*) 'Do you want to shift the grid along the x,y or z axis (y/n)?'
           write(*,*)
           shift_tau(:) = 0.0
           read (*,*) outputchar
           if (outputchar == 'y' .or. outputchar == 'Y') then
              write(*,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,') :'
              write(*,*)
              read (*,*) gridshift1, gridshift2, gridshift3
              shift_tau(:) = gridshift1*rprimd(:,1)/(nr1+1) + gridshift2*rprimd(:,2)/(nr2+1) + gridshift3*rprimd(:,3)/(nr3+1)
           end if

           allocate(filename(1))
           filename(1)=trim(output)
           write(*,*) '  The name of your data files is : '
           write(*,*) trim(filename(1)),'  for the density (norm of the wfk),'
           write(*,*)

           open(unit=unout,file=filename(1),status='replace',form='formatted')
           rewind(unout)
           do iband=1,nband(ckpt)
              write(unout,'(a,2f20.16)')'#', eigen(iband),occ1(iband)
           end do

           write(6, '(/,a,2x,3i5)' )' Number of points per side: ',nr1+1,nr2+1,nr3+1
           write(6, '(/,a,2x,i10,//)' )' Total number of points:', (nr1+1)*(nr2+1)*(nr3+1)
           write(6,*) ' znucl = ', znucl, ' typat = ', typat, ' ntypat = ', ntypat

           write(unout,'(1X,A)')  'DIM-GROUP'
           write(unout,*) '3  1'
           write(unout,'(1X,A)') 'PRIMVEC'
           do ir1 = 1,3
              write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3)
           end do
           write(unout,'(1X,A)') 'PRIMCOORD'
           write(unout,*) natom, ' 1'
           !
           !   generate translated coordinates to match density shift
           !
           do iatom = 1,natom
              tau2 (:,iatom) = tau(:,iatom) - shift_tau(:)
           end do

           do iatom = 1,natom
              write(unout,'(i9,3(3X,ES17.10))') int(znucl(typat(iatom))), &
                   &                 Bohr_Ang*tau2(1,iatom), &
                   &                 Bohr_Ang*tau2(2,iatom), &
                   &                 Bohr_Ang*tau2(3,iatom)
           end do
           write(unout,'(1X,A)') 'ATOMS'
           do iatom = 1,natom
              write(unout,'(i9,3(3X,ES17.10))') int(znucl(typat(iatom))), &
                   &                 Bohr_Ang*tau2(1,iatom), &
                   &                 Bohr_Ang*tau2(2,iatom), &
                   &                 Bohr_Ang*tau2(3,iatom)
           end do

           !           write(unout,'(1X,A)') 'FRAMES'
           write(unout,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D'
           write(unout,*) 'datagrids'
           write(unout,'(1X,A)') 'DATAGRID_3D_DENSITY'
           write(unout,*) nr1+1,nr2+1,nr3+1
           write(unout,*) '0.0 0.0 0.0 '
           do ir1 = 1,3
              write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3)
           end do

           do ir3=gridshift3+1,nr3+1
              ii3=mod(ir3-1,nr3) + 1
              do ir2=gridshift2+1,nr2+1
                 ii2=mod(ir2-1,nr2) + 1
                 do ir1=gridshift1+1,nr1+1
                    ii1=mod(ir1-1,nr1) + 1
                    tmpr=fofr(1,ii1,ii2,ii3)
                    tmpi=fofr(2,ii1,ii2,ii3)
                    write(unout,'(e10.5)') tmpr*tmpr + tmpi*tmpi
                 end do
                 do ir1=1,gridshift1
                    ii1=mod(ir1-1,nr1) + 1
                    tmpr=fofr(1,ii1,ii2,ii3)
                    tmpi=fofr(2,ii1,ii2,ii3)
                    write(unout,'(e10.5)') tmpr*tmpr + tmpi*tmpi
                 end do
              end do
              do ir2=1,gridshift2
                 ii2=mod(ir2-1,nr2) + 1
                 do ir1=gridshift1+1,nr1+1
                    ii1=mod(ir1-1,nr1) + 1
                    tmpr=fofr(1,ii1,ii2,ii3)
                    tmpi=fofr(2,ii1,ii2,ii3)
                    write(unout,'(e10.5)') tmpr*tmpr + tmpi*tmpi
                 end do
                 do ir1=1,gridshift1
                    ii1=mod(ir1-1,nr1) + 1
                    tmpr=fofr(1,ii1,ii2,ii3)
                    tmpi=fofr(2,ii1,ii2,ii3)
                    write(unout,'(e10.5)') tmpr*tmpr + tmpi*tmpi
                 end do
              end do
           end do
           do ir3=1,gridshift3
              ii3=mod(ir3-1,nr3) + 1
              do ir2=gridshift2+1,nr2+1
                 ii2=mod(ir2-1,nr2) + 1
                 do ir1=gridshift1+1,nr1+1
                    ii1=mod(ir1-1,nr1) + 1
                    tmpr=fofr(1,ii1,ii2,ii3)
                    tmpi=fofr(2,ii1,ii2,ii3)
                    write(unout,'(e10.5)') tmpr*tmpr + tmpi*tmpi
                 end do
                 do ir1=1,gridshift1
                    ii1=mod(ir1-1,nr1) + 1
                    tmpr=fofr(1,ii1,ii2,ii3)
                    tmpi=fofr(2,ii1,ii2,ii3)
                    write(unout,'(e10.5)') tmpr*tmpr + tmpi*tmpi
                 end do
              end do
              do ir2=1,gridshift2
                 ii2=mod(ir2-1,nr2) + 1
                 do ir1=gridshift1+1,nr1+1
                    ii1=mod(ir1-1,nr1) + 1
                    tmpr=fofr(1,ii1,ii2,ii3)
                    tmpi=fofr(2,ii1,ii2,ii3)
                    write(unout,'(e10.5)') tmpr*tmpr + tmpi*tmpi
                 end do
                 do ir1=1,gridshift1
                    ii1=mod(ir1-1,nr1) + 1
                    tmpr=fofr(1,ii1,ii2,ii3)
                    tmpi=fofr(2,ii1,ii2,ii3)
                    write(unout,'(e10.5)') tmpr*tmpr + tmpi*tmpi
                 end do
              end do
           end do


           write(unout,*)
           write(unout,'(1X,A)') 'END_DATAGRID_3D'
           write(unout,'(1X,A)') 'END_BLOCK_DATAGRID3D'
           close(unout)

           close(unout)
           deallocate(filename)

           write(*,*)
           exit


        case(12)
# if defined NETCDF
           allocate(partwf(nr1,nr2,nr3))
           allocate(filename(1))
           filename(1)=trim(output)//'.nc'
           write(*,*) '  The name of your data files is : '   
           write(*,*) trim(filename(1)),' is your NetCDF file' 
           write(*,*)

           !Creating NetCDF file
           status = nf90_create(filename(1), nf90_clobber, ncid)
           if (status /= nf90_noerr) call handle_ncerr(status, "Creating file")

           !Ask for a title

           write(*,*) 'Do you want a title in your NetCDF file? (0 = NO, 1 = YES)'

           read(*,*) titlechoice 
           do 
              if (titlechoice ==0 .or. titlechoice ==1) exit
              write(*,*) 'The answer is not correct, you must enter a integer between 0 and 1 (0 = NO, 1 = YES)'
              read(*,*) titlechoice
           end do
           if (titlechoice ==1) then
              write(*,*) 'Enter your file''s title'
              read(*,'(A)') filetitle
           else
              write(*,*) 'No title will be add in your NetCDF file'
           end if

           kptvar = kpt(1:3,ckpt)
           originatt(1:3,1:3)=0

           do igrid = 1,3
              gridwavefun1(igrid,1)=0
              gridwavefun2(igrid,1)=0
              gridwavefun3(igrid,1)=0
              gridwavefun1(igrid,2)=Bohr_Ang*rprimd(igrid,3)/nr3
              gridwavefun2(igrid,2)=Bohr_Ang*rprimd(igrid,2)/nr2
              gridwavefun3(igrid,2)=Bohr_Ang*rprimd(igrid,1)/nr1
           end do

           !Defining dimensions

           status = nf90_def_dim(ncid,"gridsize1",nr1, gridsize1DimID)
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining dimensions")

           status = nf90_def_dim(ncid,"gridsize2",nr2, gridsize2DimID)
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining dimensions")

           status = nf90_def_dim(ncid,"gridsize3",nr3, gridsize3DimID)
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining dimensions")

           status = nf90_def_dim(ncid, "lat",3, latDimID)
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining dimensions")

           status = nf90_def_dim(ncid, "pos",2, posDimID)
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining dimensions")

           status = nf90_def_dim(ncid, "nbatom",natom, nbatomDimID)
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining dimensions")

           !Defining variables

           status = nf90_def_var(ncid, "kpoint",nf90_float, latDimID, kpointVarID)
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining variables")

           status = nf90_def_var(ncid, "latticevec",nf90_float, (/ latDimID, latDimID /), latticevecVarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")

           status = nf90_def_var(ncid, "origin",nf90_float, (/ latDimID, latDimID /), originVarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")

           status = nf90_def_var(ncid, "atomposi",nf90_float, (/ latDimID, nbatomDimID /), atomposiVarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")

           status = nf90_def_var(ncid, "atomicnum",nf90_float, nbatomDimID, atomicnumVarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")

           status = nf90_def_var(ncid, "grid1",nf90_float, (/latDimID,posDimID/), grid1VarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")

           status = nf90_def_var(ncid, "grid2",nf90_float, (/latDimID,posDimID/), grid2VarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")

           status = nf90_def_var(ncid, "grid3",nf90_float, (/latDimID,posDimID/), grid3VarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")        
           
           status = nf90_def_var(ncid, "imagwavefunction",nf90_float, (/ gridsize1DimID, gridsize2DimID, gridsize3DimID /), imagwavefunVarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")

           status = nf90_def_var(ncid, "realwavefunction",nf90_float, (/ gridsize1DimID, gridsize2DimID, gridsize3DimID /), realwavefunVarID)
           if (status /= nf90_noerr) call  handle_ncerr(status, "Defining variables")

           !Defining attributes

           status = nf90_put_att(ncid,latticevecVarID , "field", "latticevec,vector")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,originVarID , "field", "origin,vector")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,atomposiVarID , "field","atomposi,vector")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,latticevecVarID , "positions","origin" )
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,atomicnumVarID , "positions", "atomposi")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,realwavefunVarID , "positions", "grid1,product,compact;grid2,product,compact;grid3,product,compact")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,realwavefunVarID , "connections", (/nr3,nr2,nr1/))
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,imagwavefunVarID , "positions", "grid1,product,compact;grid2,product,compact;grid3,product,compact")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,imagwavefunVarID , "connections", (/nr3,nr2,nr1/))
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,kpointVarID , "long_name", "K-point Value")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,latticevecVarID , "long_name", "Lattice Vectors")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,atomposiVarID , "long_name", "Atomic Positions")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid, atomicnumVarID, "long_name", "Atomic Numbers")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid, imagwavefunVarID, "long_name", " Imaginary Part Wave Function")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid, realwavefunVarID, "long_name", " Real Part Wave Function")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid, originVarID, "long_name", "Origin of the Lattice")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,latticevecVarID, "units", "Angstroms")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           status = nf90_put_att(ncid,atomposiVarID, "units", "Angstroms")
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           if (titlechoice ==1) then
              status = nf90_put_att(ncid, nf90_global, "title", filetitle)
              if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")
           end if

           !Add the creation date
           call date_and_time(strdat,strtime,strzone,values)
           yyyy=values(1)
           mm=values(2)
           dd=values(3)
           write(stridate(1:2),'(I2)') dd
           stridate(3:3)=" "
           stridate(4:7)=monnam(mm)
           write(stridate(8:11),'(I4)') yyyy

           status = nf90_put_att(ncid, nf90_global,"date", stridate)
           if (status /= nf90_noerr) call handle_ncerr(status, "Defining attributes")

           !Ending the define mode and entering data mode 

           status = nf90_enddef(ncid)
           if (status /= nf90_noerr) call handle_ncerr(status, "Entering data mode")

           !Defining the variables

           status = nf90_put_var(ncid,kpointVarID,kptvar)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           status = nf90_put_var(ncid,latticevecVarID,Bohr_Ang*rprimd)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           status = nf90_put_var(ncid,originVarID,originatt)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           status = nf90_put_var(ncid,atomposiVarID, Bohr_Ang*tau)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           status = nf90_put_var(ncid,atomicnumVarID,znucl(typat))
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           status = nf90_put_var(ncid,grid1VarID,gridwavefun1)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           status = nf90_put_var(ncid,grid2VarID,gridwavefun2)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           status = nf90_put_var(ncid,grid3VarID,gridwavefun3)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           partwf = fofr(2,1:nr1,1:nr2,1:nr3)

           status = nf90_put_var(ncid,imagwavefunVarID, partwf)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           partwf = fofr(1,1:nr1,1:nr2,1:nr3)

           status = nf90_put_var(ncid,realwavefunVarID, partwf)
           if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")

           status = nf90_close(ncid)
           if (status /= nf90_noerr) call handle_ncerr(status, "Closing file")


           write(*,*) 'The NetCDF file is done'
           !The NetCDF file is done

           deallocate(partwf)
           deallocate(filename)

           write(*,*)
           exit

# else 
           write(*,*) 'NetCDF is not defined. You must choose another option'
           exit
# endif

           !************************************************************

        case(13)
           write(*,*)
           write(*,*) 'Give 1 files of formatted data'
           write(*,*) 'The files are ready to be used with XCrysDen'
           write(*,*)
           gridshift1 = 0
           gridshift2 = 0
           gridshift3 = 0
           write(*,*) 'Do you want to shift the grid along the x,y or z axis (y/n)?'
           write(*,*)
           shift_tau(:) = 0.0
           read (*,*) outputchar
           if (outputchar == 'y' .or. outputchar == 'Y') then
              write(*,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,') :'
              write(*,*)
              read (*,*) gridshift1, gridshift2, gridshift3
              shift_tau(:) = gridshift1*rprimd(:,1)/(nr1+1) + gridshift2*rprimd(:,2)/(nr2+1) + gridshift3*rprimd(:,3)/(nr3+1)
           end if

           allocate(filename(1))
           filename(1)=trim(output)
           write(*,*) '  The name of your data files is : '
           write(*,*) trim(filename(1)),'  for the wavefunction,'
           write(*,*)

           open(unit=unout,file=filename(1),status='replace',form='formatted')
           rewind(unout)
           do iband=1,nband(ckpt)
              write(unout,'(a,2f20.16)')'#', eigen(iband),occ1(iband)
           end do

           write(6, '(/,a,2x,3i5)' )' Number of points per side: ',nr1+1,nr2+1,nr3+1
           write(6, '(/,a,2x,i10,//)' )' Total number of points:', (nr1+1)*(nr2+1)*(nr3+1)
           write(6,*) ' znucl = ', znucl, ' typat = ', typat, ' ntypat = ', ntypat

           write(unout,'(1X,A)')  'DIM-GROUP'
           write(unout,*) '3  1'
           write(unout,'(1X,A)') 'PRIMVEC'
           do ir1 = 1,3
              write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3)
           end do
           write(unout,'(1X,A)') 'PRIMCOORD'
           write(unout,*) natom, ' 1'
           !
           !   generate translated coordinates to match density shift
           !
           do iatom = 1,natom
              tau2 (:,iatom) = tau(:,iatom) - shift_tau(:)
           end do

           do iatom = 1,natom
              write(unout,'(i9,3(3X,ES17.10))') int(znucl(typat(iatom))), &
                   &                 Bohr_Ang*tau2(1,iatom), &
                   &                 Bohr_Ang*tau2(2,iatom), &
                   &                 Bohr_Ang*tau2(3,iatom)
           end do
           write(unout,'(1X,A)') 'ATOMS'
           do iatom = 1,natom
              write(unout,'(i9,3(3X,ES17.10))') int(znucl(typat(iatom))), &
                   &                 Bohr_Ang*tau2(1,iatom), &
                   &                 Bohr_Ang*tau2(2,iatom), &
                   &                 Bohr_Ang*tau2(3,iatom)
           end do

           !           write(unout,'(1X,A)') 'FRAMES'
           
           ! Begin the wavefunction block data
           write(unout,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D'
           write(unout,*) 'Wavefunction_Data'
           
           ! Write Real part of wavefunction
           write(unout,'(1X,A)') 'BEGIN_DATAGRID_3D_Real_Part_WFK'
           write(unout,*) nr1,nr2,nr3
           write(unout,*) '0.0 0.0 0.0 '
           do ir1 = 1,3
              write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3)
           end do
 
           !open(unit=unout,file=output,status='replace',form='formatted')
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    write(unout,'(ES17.10)') fofr(1,ir1,ir2,ir3)   
                 end do
              end do
           end do
           write(unout,*)
           write(unout,'(1X,A)') 'END_DATAGRID_3D'
           !close(unout)
           

           ! Write Imaginary part of wavefunction
           write(unout,'(1X,A)') 'BEGIN_DATAGRID_3D_Imaginary_Part_WFK'
           write(unout,*) nr1,nr2,nr3
           write(unout,*) '0.0 0.0 0.0 '
           do ir1 = 1,3
              write(unout,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(ir2,ir1), ir2=1,3)
           end do

           !open(unit=unout,file=output,status='replace',form='formatted')
           do ir3=1,nr3
              do ir2=1,nr2
                 do ir1=1,nr1
                    write(unout,'(ES17.10)') fofr(2,ir1,ir2,ir3)  
                 end do
              end do
           end do
           write(unout,*)
           write(unout,'(1X,A)') 'END_DATAGRID_3D'

           ! End wavefunction block data
           write(unout,'(1X,A)') 'END_BLOCK_DATAGRID3D'
           
           !close(unout)

           close(unout)
           deallocate(filename)

           write(*,*)
           exit


        case(14)
           write(6,*)' Exit inner loop'
           !           stop
           select_exit = 1


        case default
           write(6,*) ' This choice is not valid.'
           write(6,*)
           cycle

        end select

     end do

     ckpt=oldckpt
     cband=oldcband
     csppol=oldcsppol
     cspinor=oldcspinor
     !deallocate the datas
     deallocate(fofr)

     write(*,*) ' Task ',ichoice,' has been done !'
     write(*,*)
     write(*,*) ' Run interpolation again? (1=default=yes,2=no)'
     read(*,*) iprompt
     if(iprompt==2) then
        exit
     else
        cycle
     end if
  end do
  !deallocate the datas
  deallocate(cg,eigen,kg_dum,ph1d,occ1)

  ! Close the WF file
  close(19)

end subroutine wffile
!!***
