TABLE OF CONTENTS


ABINIT/m_iowf [ Modules ]

[ Top ] [ Modules ]

NAME

 m_iowf

FUNCTION

 Procedures for the IO of the WFK file.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, AR, MB, MVer, MG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_iowf
24 
25  use defs_basis
26  use defs_abitypes
27  use defs_wvltypes
28  use m_abicore
29  use m_errors
30  use m_xmpi
31  use m_wffile
32  use m_abi_etsf
33  use m_nctk
34  use m_wfk
35 #ifdef HAVE_NETCDF
36  use netcdf
37 #endif
38  use m_hdr
39  use m_ebands
40 
41  use m_time,           only : cwtime, timab
42  use m_io_tools,       only : get_unit, flush_unit, iomode2str
43  use m_fstrings,       only : endswith, sjoin
44  use m_numeric_tools,  only : mask2blocks
45  use defs_datatypes,   only : ebands_t, pseudopotential_type
46  use m_cgtools,        only : cg_zcopy
47  use m_crystal,        only : crystal_t, crystal_free
48  use m_crystal_io,     only : crystal_ncwrite, crystal_from_hdr
49  use m_rwwf,           only : rwwf
50  use m_mpinfo,         only : proc_distrb_cycle
51  use m_vkbr,           only : calc_vkb
52  use m_wvl_rwwf,       only : wvl_write
53 
54  implicit none
55 
56  private
57 
58  public :: outwf

m_iowf/cg2seqblocks [ Functions ]

[ Top ] [ m_iowf ] [ Functions ]

NAME

 cg2seqblocks

FUNCTION

INPUTS

  npwtot_k
  npw_k=
  cg_k(2,npw_k*nband_k)
  gmpi2seq
  comm_fft=FFT communicator

OUTPUT

  bstart
  bcount

PARENTS

      m_iowf

CHILDREN

      xmpi_sum_master

SOURCE

1578 subroutine cg2seqblocks(npwtot_k,npw_k,nband,cg_k,gmpi2seq,comm_bandfft,bstart,bcount,my_cgblock)
1579 
1580 
1581 !This section has been created automatically by the script Abilint (TD).
1582 !Do not modify the following lines by hand.
1583 #undef ABI_FUNC
1584 #define ABI_FUNC 'cg2seqblocks'
1585 !End of the abilint section
1586 
1587  implicit none
1588 
1589 !Arguments ------------------------------------
1590 !scalars
1591  integer,intent(in) :: npwtot_k,npw_k,nband,comm_bandfft
1592  integer,intent(out) :: bstart,bcount
1593 !arrays
1594  integer,intent(in) :: gmpi2seq(npw_k)
1595  real(dp),intent(in) :: cg_k(2,npw_k,nband)
1596  real(dp),allocatable,intent(out) :: my_cgblock(:,:,:)
1597 
1598 !Local variables-------------------------------
1599 !scalars
1600  integer :: me,nprocs,rank,ig,ib,igseq,ierr,nbb,nbmax,band
1601 !arrays
1602  integer,allocatable :: bstart_rank(:)
1603  real(dp),allocatable :: cgbuf(:,:,:)
1604 
1605 ! *************************************************************************
1606 
1607  me = xmpi_comm_rank(comm_bandfft); nprocs = xmpi_comm_size(comm_bandfft)
1608 
1609  ! Handle sequential case.
1610  if (nprocs == 1) then
1611    bstart = 1; bcount = nband
1612    ABI_STAT_MALLOC(my_cgblock, (2, npwtot_k, nband), ierr)
1613    ABI_CHECK(ierr==0, "oom my_cgblock")
1614    my_cgblock = cg_k
1615    return ! DOH
1616  end if
1617 
1618  ABI_MALLOC(bstart_rank, (0:nprocs))
1619  bstart_rank = [(1 + (nband/nprocs)*rank, rank=0,nprocs-1), 1 + nband]
1620 
1621  ! Allocate MPI buffer (same size on each MPI proc)
1622  nbmax = 0
1623  do rank=0,nprocs-1
1624    nbmax = max(nbmax, bstart_rank(rank+1) - bstart_rank(rank))
1625  end do
1626  ABI_STAT_MALLOC(cgbuf, (2, npwtot_k, nbmax), ierr)
1627  ABI_CHECK(ierr==0, "oom cgbuf")
1628 
1629  nbb = bstart_rank(me+1) - bstart_rank(me)
1630  ABI_STAT_MALLOC(my_cgblock, (2, npwtot_k, nbb), ierr)
1631  ABI_CHECK(ierr==0, "oom my_cgblock")
1632 
1633  ! TODO: This should be replaced by gatherv but premature optimization....
1634  do rank=0,nprocs-1
1635    bstart = bstart_rank(rank)
1636    bcount = bstart_rank(rank+1) - bstart_rank(rank)
1637 
1638    do band=bstart, bstart+bcount-1
1639      ib = band - bstart + 1
1640      cgbuf(:,:,ib) = zero
1641      do ig=1,npw_k
1642        igseq = gmpi2seq(ig)
1643        cgbuf(:,igseq,ib) = cg_k(:,ig,band)
1644      end do
1645    end do ! band
1646 
1647    call xmpi_sum_master(cgbuf,rank,comm_bandfft,ierr)
1648    if (me == rank) my_cgblock = cgbuf(:,:,:bcount)
1649  end do ! rank
1650 
1651  bstart = bstart_rank(me)
1652  bcount = bstart_rank(me+1) - bstart_rank(me)
1653 
1654  ABI_FREE(cgbuf)
1655  ABI_FREE(bstart_rank)
1656 
1657 end subroutine cg2seqblocks

m_iowf/cg_ncwrite [ Functions ]

[ Top ] [ m_iowf ] [ Functions ]

NAME

 cg_ncwrite

FUNCTION

 Conduct output of a "wave-functions" file with netcdf

INPUTS

  fname= character string giving the root to form the name of the
   output WFK or WFQ file if response==0, otherwise it is the filename.
  dtset <type(dataset_type)>=all input variables for this dataset
  hdr <type(hdr_type)>=the header of wf, den and pot files
  response: if == 0, GS wavefunctions , if == 1, RF wavefunctions
  mpw=maximum number of plane waves
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mkmem=maximum number of k-points treated by this node
  nkpt=number of k points
  nsppol=1 for unpolarized, 2 for spin-polarized
  eigen((2*mband)**response *mband*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point
  occ(mband*nkpt*nsppol)=occupations for all bands at each k point
  cg(2,mcg)=wavefunction array (storage if nkpt>1)
  npwarr(nkpt)=number of planewaves in basis and on boundary for each k
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mpi_enreg=information about MPI parallelization

OUTPUT

  done=.True if cg_ncwrite can handle the output of the WFK file in parallel.

PARENTS

      m_iowf

CHILDREN

      xmpi_sum_master

SOURCE

 738 #ifdef HAVE_NETCDF
 739 
 740 subroutine cg_ncwrite(fname,hdr,dtset,response,mpw,mband,nband,nkpt,nsppol,nspinor,mcg,&
 741                       mkmem,eigen,occ,cg,npwarr,kg,mpi_enreg,done)
 742 
 743 
 744 !This section has been created automatically by the script Abilint (TD).
 745 !Do not modify the following lines by hand.
 746 #undef ABI_FUNC
 747 #define ABI_FUNC 'cg_ncwrite'
 748 !End of the abilint section
 749 
 750  implicit none
 751 
 752 !Arguments ------------------------------------
 753 !scalars
 754  integer,intent(in) :: response,mband,mcg,mkmem,mpw,nkpt,nsppol,nspinor
 755  character(len=*),intent(in) :: fname
 756  logical,intent(out) :: done
 757  type(dataset_type),intent(in) :: dtset
 758  type(MPI_type),intent(in) :: mpi_enreg
 759  type(hdr_type),intent(in) :: hdr
 760 !arrays
 761  integer, intent(in) :: nband(nkpt*nsppol),kg(3,mpw*mkmem),npwarr(nkpt)
 762  real(dp),intent(in) :: cg(2,mcg)
 763  real(dp),intent(in) :: eigen((2*mband)**response*mband*nkpt*nsppol),occ(mband*nkpt*nsppol)
 764 
 765 !Local variables-------------------------------
 766 !scalars
 767  integer,parameter :: master=0,fform2=2
 768  integer :: ii,iomode,icg,iband,ikg,ikpt,spin,me_cell,me_kpt,me_band,me_spinor,my_nspinor,nband_k,npw_k
 769  integer :: comm_cell,comm_fft,comm_bandfft,formeig
 770  integer :: cnt,min_cnt,max_cnt,ierr,action,source,ncid,ncerr,cg_varid,kg_varid !,eig_varid,
 771  integer :: timrev,paral_kgb,npwtot_k !,start_pwblock !,start_cgblock !count_pwblock,
 772  integer :: ipw,ispinor_index,npwso,npwsotot,npwtot,nspinortot,ikpt_this_proc,ispinor
 773  integer :: bandpp,nproc_band,nproc_fft,nproc_spinor,me_fft,nproc_cell,nwrites
 774  integer :: comm_mpiio,nranks,bstart,bcount !nbdblock,nblocks,
 775  !integer :: band_blocksize,band
 776  real(dp) :: cpu,wall,gflops
 777  logical :: ihave_data,iam_master,single_writer,same_layout,use_collective
 778  character(len=500) :: msg
 779  character(len=fnlen) :: path
 780  type(wfk_t) :: wfk
 781  type(crystal_t) :: crystal
 782  type(ebands_t) :: gs_ebands
 783 !arrays
 784  integer,allocatable :: kg_k(:,:),iter2kscgkg(:,:),ind_cg_mpi_to_seq(:),rank_has_cg(:),ranks_io(:)!,gblock(:,:)
 785  real(dp) :: tsec(2)
 786  real(dp),allocatable :: eigen3d(:,:,:),occ3d(:,:,:),cg_k(:,:),my_cgblock(:,:,:)
 787 
 788 ! *************************************************************************
 789 
 790  DBG_ENTER("COLL")
 791  done = .False.
 792 
 793  path = nctk_ncify(fname)
 794  call wrtout(std_out, sjoin("in cg_ncwrite with path:", path), "COLL")
 795 
 796  ! communicators and ranks
 797  comm_cell = mpi_enreg%comm_cell; me_cell = mpi_enreg%me_cell; nproc_cell = mpi_enreg%nproc_cell
 798  comm_fft = mpi_enreg%comm_fft; me_fft  = mpi_enreg%me_fft; nproc_fft  = mpi_enreg%nproc_fft
 799  comm_bandfft = mpi_enreg%comm_bandfft
 800  me_kpt = mpi_enreg%me_kpt; me_band = mpi_enreg%me_band; me_spinor = mpi_enreg%me_spinor
 801  iam_master = (me_kpt == master)
 802 
 803  paral_kgb = dtset%paral_kgb
 804  nproc_band = mpi_enreg%nproc_band
 805  bandpp     = mpi_enreg%bandpp
 806  nproc_spinor = mpi_enreg%nproc_spinor
 807 
 808  ! FIXME
 809  my_nspinor = max(1, nspinor/nproc_spinor)
 810  if (nspinor == 2 .and. my_nspinor == 1) then
 811    MSG_ERROR("Spinor parallelization not coded yet")
 812  end if
 813 
 814  ! TODO: Be careful with response == 1 in parallel because the distribution of the cg
 815  ! can be **VERY** different from cg if nprocs > nkpt * nsppol
 816  !ABI_CHECK(response==0, "response == 1 not coded")
 817  if (size(hdr%nband) == size(nband)) then
 818    ABI_CHECK(all(Hdr%nband == nband),"nband")
 819  else
 820    MSG_ERROR("hdr%nband and nband have different size!")
 821  end if
 822 
 823  if (xmpi_comm_size(comm_cell) == 1) then
 824    ABI_CHECK(all(npwarr == hdr%npwarr), "npwarr != hdr%npwarr")
 825  end if
 826 
 827  ! FIXME: Use abinit convention for timrev
 828  timrev = 2
 829  call crystal_from_hdr(crystal, hdr, timrev)
 830 
 831  ! TODO
 832  ! Be careful with response == 1.
 833  ! gs_ebands contains the GS eigenvalues and occupation and will be written if this is a
 834  ! GS wfk. If we have a DFPT file, eigen stores the GKK matrix element, in this case
 835  ! we don't write gs_ebands but we define new variables in the netcdf file to store the GKK
 836  ABI_MALLOC(occ3d, (mband,nkpt,nsppol))
 837  call unpack_eneocc(nkpt,nsppol,mband,nband,occ,occ3d)
 838 
 839  if (response == 0) then
 840    formeig = 0
 841    ABI_MALLOC(eigen3d, (mband,nkpt,nsppol))
 842    call unpack_eneocc(nkpt,nsppol,mband,nband,eigen,eigen3d)
 843    gs_ebands = ebands_from_hdr(hdr, mband, eigen3d); gs_ebands%occ = occ3d
 844    ABI_FREE(eigen3d)
 845  else
 846    formeig = 1
 847  end if
 848 
 849  ! same_layout is set to True if the interal representation of the cgs
 850  ! is compatible with the representation on file.
 851  iomode = IO_MODE_ETSF
 852  if (response == 0) then
 853     same_layout = (paral_kgb == 0 .or. (paral_kgb == 1 .and. all([nproc_fft, nproc_band, nproc_spinor] == 1)))
 854  else
 855    ! For the time being, these cases are not implemented in the DFPT part.
 856    ABI_CHECK(nproc_fft==1, "nproc_fft != 1 not coded")
 857    ABI_CHECK(nproc_band==1, "nproc_band != 1 not coded")
 858    ABI_CHECK(nproc_spinor==1, "nproc_spinor != 1 not coded")
 859 
 860    ! Note: It would be possible to use collective IO if the cg1 are block-distributed
 861    same_layout = .True.
 862    spin_loop: do spin=1,nsppol
 863      do ikpt=1,nkpt
 864        nband_k = nband(ikpt + (spin-1)*nkpt)
 865        if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then
 866          if (any(mpi_enreg%proc_distrb(ikpt,:nband_k,spin) /= me_kpt)) then
 867            same_layout = .False.; exit spin_loop
 868          end if
 869        end if
 870      end do
 871    end do spin_loop
 872    call xmpi_land(same_layout, comm_cell)
 873  end if
 874 
 875  call cwtime(cpu,wall,gflops,"start")
 876 
 877  if (same_layout) then
 878    single_writer = .True.
 879    if (nctk_has_mpiio) single_writer = .False.
 880    !single_writer = .True.
 881 
 882    write(msg,'(5a,l1)')&
 883      "same layout Writing WFK file: ",trim(path),", with iomode: ",trim(iomode2str(iomode)),", single writer: ",single_writer
 884    call wrtout(std_out,msg,'PERS', do_flush=.True.)
 885 
 886    if (.not. single_writer) then
 887 
 888 #ifdef HAVE_NETCDF
 889      ! master opens the file and write the metadata.
 890      if (xmpi_comm_rank(comm_cell) == master) then
 891        ncerr = nf90_einval
 892 #ifdef HAVE_NETCDF_MPI
 893        ncerr = nf90_create(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), &
 894          comm=xmpi_comm_self, info=xmpio_info, ncid=ncid)
 895 #endif
 896        NCF_CHECK_MSG(ncerr, sjoin("create_par: ", path))
 897 
 898        call wfk_ncdef_dims_vars(ncid, hdr, fform2, write_hdr=.True.)
 899        NCF_CHECK(crystal_ncwrite(crystal, ncid))
 900 
 901        if (response == 0) then
 902          NCF_CHECK(ebands_ncwrite(gs_ebands, ncid))
 903        else
 904          call ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d)
 905        end if
 906 
 907        NCF_CHECK(nf90_close(ncid))
 908      end if
 909 
 910      ! Compute the table for collective IO.
 911      ABI_CALLOC(iter2kscgkg, (4, nkpt*nsppol))
 912      cnt = 0; icg = 0
 913      do spin=1,nsppol
 914        ikg = 0
 915        do ikpt=1,nkpt
 916          nband_k = nband(ikpt + (spin-1)*nkpt)
 917          npw_k = npwarr(ikpt)
 918          if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then
 919            ! FIXME v2[01], v3[6], with 4 procs
 920            ! v6[7] and v7[68], v7[69] fail but due to a extra line with nprocs
 921            ! v7[96] with np=4 seems to be more serious but it crashes also in trunk.
 922            ! v2[88] is slow likely due to outscfv
 923            ABI_CHECK(all(mpi_enreg%proc_distrb(ikpt,:nband_k,spin) == me_kpt), "bands are distributed")
 924            cnt = cnt + 1
 925            iter2kscgkg(:,cnt) = [ikpt, spin, icg, ikg]
 926            icg = icg + npw_k*my_nspinor*nband_k
 927            ikg = ikg + npw_k
 928          end if
 929        end do
 930      end do
 931      if (cnt == 0) then
 932        write(std_out,*)"cnt == 0 for me_cell, me, me_kpt",mpi_enreg%me_cell, mpi_enreg%me, mpi_enreg%me_kpt
 933        !ABI_CHECK(cnt > 0, "This processor does not have wavefunctions!")
 934      end if
 935 
 936      call xmpi_min(cnt, min_cnt, comm_cell, ierr)
 937      call xmpi_max(cnt, max_cnt, comm_cell, ierr)
 938 
 939      ! Handle idle procs, i.e. processors that do not have wavefunctions
 940      ! This happens if paral_kgb == 0 and nprocs > nkpt * nsppol (Abinit does not stop anymore!)
 941      comm_mpiio = comm_cell
 942 
 943      if (min_cnt <= 0) then
 944        MSG_COMMENT("Will create subcommunicator to exclude idle processors from MPI-IO collective calls")
 945        ABI_CHECK(paral_kgb == 0, "paral_kgb == 1 with idle processors should never happen")
 946 
 947        ! Prepare the call to xmpi_subcomm that will replace comm_mpiio.
 948        ABI_CALLOC(rank_has_cg, (0:nproc_cell-1))
 949        do spin=1,nsppol
 950          do ikpt=1,nkpt
 951            nband_k = nband(ikpt + (spin-1)*nkpt)
 952            if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then
 953               rank_has_cg(me_kpt) = 1
 954               exit
 955             end if
 956          end do
 957        end do
 958 
 959        call xmpi_sum(rank_has_cg,comm_cell,ierr)
 960        nranks = count(rank_has_cg == 1)
 961        ABI_MALLOC(ranks_io, (nranks))
 962        cnt = 0
 963        do ii=0,nproc_cell-1
 964          if (rank_has_cg(ii) == 1) then
 965            cnt = cnt + 1
 966            ranks_io(cnt) = ii
 967          end if
 968        end do
 969        !write(std_out,*)"nranks, ranks_io:", nranks, ranks_io
 970        comm_mpiio = xmpi_subcomm(comm_cell, nranks, ranks_io)
 971        if (.not. rank_has_cg(me_kpt) == 1) then
 972          comm_mpiio = xmpi_comm_null
 973          ABI_CHECK(rank_has_cg(me_kpt) == 0, "rank_has_cg must be 0 or 1")
 974        end if
 975        ABI_FREE(ranks_io)
 976        ABI_FREE(rank_has_cg)
 977      end if
 978 
 979      ! Open the file in parallel inside comm_mpiio.
 980      call xmpi_barrier(comm_cell)
 981      if (comm_mpiio == xmpi_comm_null) goto 100
 982 
 983      ncerr = nf90_einval
 984 #ifdef HAVE_NETCDF_MPI
 985      ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write),&
 986                        comm=comm_mpiio, info=xmpio_info, ncid=ncid)
 987 #endif
 988      NCF_CHECK_MSG(ncerr, sjoin("open_par: ", path))
 989 
 990      ! Use individual IO (default) for the G-vectors [3, mpw, nkpt]
 991      NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid))
 992      spin = 1; ikg=0
 993      do ikpt=1,nkpt
 994        npw_k = npwarr(ikpt)
 995        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,1,spin,me_kpt)) cycle
 996        ncerr = nf90_put_var(ncid, kg_varid, kg(:, 1+ikg:), start=[1,1,ikpt], count=[3,npw_k,1])
 997        NCF_CHECK_MSG(ncerr, "put_kg_k")
 998        ikg = ikg + npw_k
 999      end do
1000 
1001      NCF_CHECK(nf90_inq_varid(ncid, "coefficients_of_wavefunctions", cg_varid))
1002 
1003      use_collective = (response == 0) ! or (response == 1 .and. nproc_cell == 1)
1004 
1005      if (use_collective) then
1006        call wrtout(std_out,"Using collective IO for the CGs","COLL")
1007        ! Use collective IO for the CGs
1008        ncerr = nf90_einval
1009 #ifdef HAVE_NETCDF_MPI
1010        ncerr = nf90_var_par_access(ncid, cg_varid, nf90_collective)
1011 #endif
1012        NCF_CHECK(ncerr)
1013 
1014        do cnt=1,max_cnt
1015           ikpt = iter2kscgkg(1,cnt)
1016           spin = iter2kscgkg(2,cnt)
1017           icg  = iter2kscgkg(3,cnt)
1018           ikg  = iter2kscgkg(4,cnt)
1019 
1020           ! The array on file has shape [cplex, mpw, nspinor, mband, nkpt, nsppol]
1021           if (ikpt /= 0) then
1022             nband_k = nband(ikpt + (spin-1)*nkpt)
1023             npw_k   = npwarr(ikpt)
1024 
1025             !write(std_out,*)"About to write ikpt, spin, npw_k, icg: ",ikpt, spin, npw_k, icg
1026             ncerr = nf90_put_var(ncid, cg_varid, cg(:, 1+icg:), start=[1,1,1,1,ikpt,spin], &
1027                count=[2,npw_k,nspinor,nband_k,1,1])
1028             NCF_CHECK_MSG(ncerr, "writing cg")
1029           else
1030             ! This happens when nkpt * nsppol // nprocs != 0
1031             ! Note that we are using collective MPI-IO hence all processors must call put_var
1032             ! Here we re-write the ug(0) of the first (k-point, spin) treated by the node.
1033             ikpt = iter2kscgkg(1,1)
1034             spin = iter2kscgkg(2,1)
1035             ncerr = nf90_put_var(ncid, cg_varid, cg, start=[1,1,1,1,ikpt,spin], count=[1,1,1,1,1,1])
1036             NCF_CHECK_MSG(ncerr, "re-writing cg")
1037           end if
1038        end do
1039 
1040      else
1041        call wrtout(std_out,"Using individual IO for the CGs","COLL")
1042        ! Individual IO of the CGs (for debugging purposes)
1043        icg = 0
1044        do spin=1,nsppol
1045          do ikpt=1,nkpt
1046            nband_k = nband(ikpt + (spin-1)*nkpt)
1047            npw_k = npwarr(ikpt)
1048            if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) cycle
1049            ncerr = nf90_put_var(ncid, cg_varid, cg(:, 1+icg:), start=[1,1,1,1,ikpt,spin], count=[2,npw_k,nspinor,nband_k,1,1])
1050            NCF_CHECK_MSG(ncerr, "put_var")
1051            icg = icg+npw_k*my_nspinor*nband_k
1052          end do
1053        end do
1054      end if
1055 
1056      NCF_CHECK(nf90_close(ncid))
1057 
1058 100  call xmpi_barrier(comm_cell)
1059      ABI_FREE(iter2kscgkg)
1060 
1061      done = .True.
1062      call cwtime(cpu,wall,gflops,"stop")
1063      write(msg,'(2(a,f8.2))')" collective ncwrite, cpu: ",cpu,", wall: ",wall
1064      call wrtout(std_out,msg,"PERS")
1065 #endif
1066 ! HAVE_NETCDF
1067    else ! single_writer
1068      if (nproc_cell > 1) then
1069        MSG_WARNING("Slow version without MPI-IO support. Processors send data to master...")
1070      else
1071        call wrtout(std_out, "Using netcdf library without MPI-IO support")
1072      end if
1073 
1074      ABI_MALLOC(kg_k,(3,mpw))
1075      ABI_STAT_MALLOC(cg_k,(2,mpw*my_nspinor*mband), ierr)
1076      ABI_CHECK(ierr==0, "out of memory in cg_k")
1077 
1078      if (iam_master) then
1079        call wfk_open_write(wfk,hdr,path,formeig,iomode,get_unit(),xmpi_comm_self,write_hdr=.True.)
1080 
1081        NCF_CHECK(crystal_ncwrite(crystal, wfk%fh))
1082        !write(std_out,*)"after crystal_ncwrite"
1083 
1084        ! Write eigenvalues and occupations (these arrays are not MPI-distributed)
1085        if (response == 0) then
1086          NCF_CHECK(ebands_ncwrite(gs_ebands, wfk%fh))
1087        else
1088          call ncwrite_eigen1_occ(wfk%fh, nband, mband, nkpt, nsppol, eigen, occ3d)
1089        end if
1090      end if
1091 
1092      icg = 0
1093      do spin=1,nsppol
1094        ikg = 0
1095        do ikpt=1,nkpt
1096          nband_k = nband(ikpt + (spin-1)*nkpt)
1097          npw_k   = npwarr(ikpt)
1098 
1099          call xmpi_barrier(comm_cell)
1100 
1101          ! Transfer the wavefunctions and the g-vectors to the master processor
1102          source = minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,spin))
1103          ihave_data = (source==me_kpt)
1104 
1105          action=0
1106          if (iam_master .and. ihave_data)    action=1 ! I am the master node, and I have the data in cg
1107          if (.not.iam_master.and.ihave_data) action=2 ! I am not the master, and I have the data => send to master
1108          if (iam_master.and..not.ihave_data) action=3 ! I am the master, and I receive the data
1109 
1110          if (action==1) then ! Copy from kg and cg
1111            kg_k(:,1:npw_k) = kg(:,ikg+1:ikg+npw_k)
1112            call cg_zcopy(npw_k*my_nspinor*nband_k, cg(1,icg+1), cg_k)
1113          end if
1114 
1115          ! Exchange data
1116          if (action==2.or.action==3) then
1117            call timab(48,1,tsec)
1118            if (action==2) then
1119              call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,source,kg_k,master,comm_cell,ierr)
1120              call xmpi_exch(cg(:,icg+1:icg+nband_k*npw_k*my_nspinor),2*nband_k*npw_k*my_nspinor,&
1121 &             source,cg_k,master,comm_cell,ierr)
1122            else
1123              call xmpi_exch(kg_k,3*npw_k,source,kg_k,master,comm_cell,ierr)
1124              call xmpi_exch(cg_k,2*nband_k*npw_k*my_nspinor,source,cg_k,master,comm_cell,ierr)
1125            end if
1126            call timab(48,2,tsec)
1127          end if
1128 
1129          ! Master writes this block of bands.
1130          if (iam_master) then
1131            if (response == 0) then
1132              call wfk_write_band_block(wfk,[1,nband_k],ikpt,spin,xmpio_single,kg_k=kg_k,cg_k=cg_k,&
1133                eig_k=gs_ebands%eig(:,ikpt,spin),occ_k=gs_ebands%occ(:,ikpt,spin))
1134              !write(std_out,*)"cg_k",cg_k(:,1:2)
1135            else
1136              call wfk_write_band_block(wfk,[1,nband_k],ikpt,spin,xmpio_single,kg_k=kg_k,cg_k=cg_k)
1137            end if
1138          end if
1139 
1140          if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt))) then
1141            icg = icg+npw_k*my_nspinor*nband_k
1142            ikg = ikg+npw_k
1143          end if
1144 
1145        end do !ikpt
1146      end do !spin
1147 
1148      ABI_FREE(kg_k)
1149      ABI_FREE(cg_k)
1150 
1151      if (iam_master) call wfk_close(wfk)
1152      call xmpi_barrier(comm_cell)
1153 
1154      done = .True.
1155      call cwtime(cpu,wall,gflops,"stop")
1156      write(msg,'(2(a,f8.2))')" individual ncwrite, cpu: ",cpu,", wall: ",wall
1157      call wrtout(std_out,msg,"PERS")
1158    end if
1159 
1160  else ! not same_layout
1161 
1162    if (nctk_has_mpiio) then
1163 #ifdef HAVE_NETCDF
1164      call wrtout(std_out, &
1165        sjoin("scattered data. writing WFK file",trim(path),", with iomode",iomode2str(iomode)), 'PERS', do_flush=.True.)
1166 
1167      ! master write the metadata.
1168      if (xmpi_comm_rank(comm_cell) == master) then
1169        ncerr = nf90_einval
1170 #ifdef HAVE_NETCDF_MPI
1171        ncerr = nf90_create(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), &
1172          comm=xmpi_comm_self, info=xmpio_info, ncid=ncid)
1173 #endif
1174        NCF_CHECK_MSG(ncerr, sjoin("create_par:", path))
1175 
1176        call wfk_ncdef_dims_vars(ncid, hdr, fform2, write_hdr=.True.)
1177        NCF_CHECK(crystal_ncwrite(crystal, ncid))
1178 
1179        ! Write eigenvalues and occupations (these arrays are not MPI-distributed)
1180        if (response == 0) then
1181          NCF_CHECK(ebands_ncwrite(gs_ebands, ncid))
1182        else
1183          call ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d)
1184        end if
1185 
1186        NCF_CHECK(nf90_close(ncid))
1187      end if
1188 
1189      ! Reopen the file inside comm_cell
1190      call xmpi_barrier(comm_cell)
1191      ncerr = nf90_einval
1192 #ifdef HAVE_NETCDF_MPI
1193      ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), &
1194        comm=comm_cell, info=xmpio_info, ncid=ncid)
1195 #endif
1196      NCF_CHECK_MSG(ncerr, sjoin("create_par:", path))
1197 
1198      ! Get var ids
1199      NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid))
1200      NCF_CHECK(nf90_inq_varid(ncid, "coefficients_of_wavefunctions", cg_varid))
1201 
1202      ! Write the G-vectors
1203      ikg = 0
1204      do ikpt=1,nkpt
1205        npw_k = npwarr(ikpt)
1206        npwtot_k = hdr%npwarr(ikpt)
1207        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,1,1,me_kpt)) cycle
1208        if (me_spinor /= 0) cycle
1209        !write(std_out,*)"In G-vector loop ",ikpt,", with me_cell me_kpt me_band, me_spinor ",me_cell,me_kpt,me_band,me_spinor
1210        !write(std_out,*)" with ik, npw_k, npwtot_k: ",ikpt, npw_k, npwtot_k
1211 
1212        ABI_MALLOC(ind_cg_mpi_to_seq, (npw_k))
1213        if (allocated(mpi_enreg%my_kgtab)) then
1214          ikpt_this_proc = mpi_enreg%my_kpttab(ikpt)
1215          ind_cg_mpi_to_seq = mpi_enreg%my_kgtab(1:npw_k,ikpt_this_proc)
1216        else
1217          ABI_CHECK(nproc_fft==1, "nproc_fft !=1 and my_kgtab not allocated")
1218          ind_cg_mpi_to_seq(1:npw_k) = [(ipw, ipw=1,npw_k)]
1219        end if
1220 
1221        ABI_CALLOC(kg_k, (3, npwtot_k))
1222        do ipw=1,npw_k
1223          kg_k(:, ind_cg_mpi_to_seq(ipw)) = kg(:,ikg+ipw)
1224        end do
1225        call xmpi_sum_master(kg_k,master,comm_bandfft,ierr)
1226        if (xmpi_comm_rank(comm_bandfft) == master) then
1227          ncerr = nf90_put_var(ncid, kg_varid, kg_k, start=[1,1,ikpt], count=[3,npwtot_k,1])
1228          NCF_CHECK_MSG(ncerr, "putting kg_k")
1229        end if
1230        ABI_FREE(kg_k)
1231 
1232        ! gblock contains block-distributed G-vectors inside comm_fft
1233        !call kg2seqblocks(npwtot_k,npw_k,kg(:,ikg+1:),ind_cg_mpi_to_seq,comm_fft,start_pwblock,count_pwblock,gblock)
1234        !write(std_out,*)"gblock(:,2)",gblock(:,2)
1235        !ncerr = nf90_put_var(ncid, kg_varid, gblock, start=[1,start_pwblock,ikpt], count=[3,count_pwblock,1])
1236        !NCF_CHECK_MSG(ncerr, "putting kg_k")
1237        !ABI_FREE(gblock)
1238 
1239        ABI_FREE(ind_cg_mpi_to_seq)
1240 
1241        ikg = ikg+npw_k
1242      end do
1243 
1244      ! Write wavefunctions
1245      if (response == 1) then
1246        ! The cg1 is allocated with size mcg1=mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol (see dfpt_looppert)
1247        ! hence bands are not MPI distributed, this is the reason why we have to use the loop over bands
1248        ! and the check on nwrites
1249        icg = 0
1250        do spin=1,nsppol
1251          do ikpt=1,nkpt
1252            nband_k = nband(ikpt + (spin-1)*nkpt)
1253            npw_k = npwarr(ikpt)
1254            npwtot_k = hdr%npwarr(ikpt)
1255            nwrites = 0
1256            do iband=1,nband_k
1257               if (mpi_enreg%proc_distrb(ikpt,iband,spin)/=me_kpt) cycle
1258               ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol]
1259               nwrites = nwrites + 1
1260               ii = 1 + (iband-1)*npw_k*my_nspinor + icg
1261               ncerr = nf90_put_var(ncid, cg_varid, cg(1:,ii:), start=[1,1,1,iband,ikpt,spin], count=[2,npwtot_k,nspinor,1,1,1])
1262               NCF_CHECK_MSG(ncerr, "put_var cg")
1263            end do ! iband
1264            if (nwrites /= 0) icg = icg + npw_k*my_nspinor*nband_k
1265          end do !ikpt
1266        end do !spin
1267 
1268      else
1269        icg = 0
1270        do spin=1,nsppol
1271          do ikpt=1,nkpt
1272            nband_k = nband(ikpt + (spin-1)*nkpt)
1273            npw_k = npwarr(ikpt)
1274            npwtot_k = hdr%npwarr(ikpt)
1275 
1276            if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) cycle
1277            !write(std_out,*)"In u(g)-vector loop ",ikpt,", with me_cell me_kpt me_band, me_spinor ",me_cell,me_kpt,me_band,me_spinor
1278            !write(std_out,*)"nband_k, npw_k, npwtot_k: ",nband_k, npw_k, npwtot_k
1279 
1280            ! Part taken from writewf
1281            ! Note that ind_cg_mpi_to_seq is wrong in nspinor > 1
1282            npwtot=npw_k; npwso=npw_k*nspinor
1283            npwsotot=npwso
1284            nspinortot = min(2,(1+mpi_enreg%paral_spinor)*nspinor)
1285 
1286            ABI_MALLOC(ind_cg_mpi_to_seq, (npwso))
1287            if (allocated(mpi_enreg%my_kgtab)) then
1288              ikpt_this_proc=mpi_enreg%my_kpttab(ikpt)
1289              do ispinor=1,nspinor
1290                ispinor_index=ispinor
1291                if (nproc_spinor > 1) ispinor_index = mpi_enreg%me_spinor + 1
1292                ind_cg_mpi_to_seq(1+npw_k*(ispinor-1):npw_k*ispinor)=npwtot*(ispinor_index-1) &
1293 &               + mpi_enreg%my_kgtab(1:npw_k,ikpt_this_proc)
1294              end do
1295            else
1296              ABI_CHECK(nproc_fft==1, "nproc_fft !=1 and my_kgtab not allocated")
1297              ind_cg_mpi_to_seq(1:npwso) = [(ipw, ipw=1,npwso)]
1298            end if
1299 
1300            ! TODO: Blocking + collective IO in comm_cell
1301            !band_blocksize = nband_k/nbdblock
1302            !ibandmin = 1
1303            !step=min(ii,MAXBAND, nband_disk)
1304            !do iblock=1,nband_disk/step+1
1305            !   ibandmax = min(ibandmin+step-1, nband_disk)
1306            !   nband_block = ibandmax-ibandmin+1
1307            !gbase = icg + npw_k*my_nspinor*nband_k
1308 
1309            ! Each MPI proc write bcount bands with npwtot_k G-vectors starting from bstart.
1310            call cg2seqblocks(npwtot_k,npw_k,nband_k,cg(:,icg+1:),ind_cg_mpi_to_seq,comm_bandfft,bstart,bcount,my_cgblock)
1311            !write(std_out,*)"bstart, bcount",bstart,bcount
1312            !write(std_out,*)"cg(:,ig+1)",cg(:,icg+1)
1313            !write(std_out,*)"my_cgblock(:,ig+1)",my_cgblock(:,1,1)
1314 
1315            ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol]
1316            ncerr = nf90_put_var(ncid, cg_varid, my_cgblock, start=[1,1,1,bstart,ikpt,spin], count=[2,npwtot_k,nspinor,bcount,1,1])
1317            NCF_CHECK_MSG(ncerr, "put_var cg")
1318 
1319            ABI_FREE(my_cgblock)
1320            ABI_FREE(ind_cg_mpi_to_seq)
1321 
1322            icg = icg+npw_k*my_nspinor*nband_k
1323          end do !ikpt
1324        end do !spin
1325      end if
1326 
1327      NCF_CHECK(nf90_close(ncid))
1328      done = .True.
1329 
1330      call cwtime(cpu,wall,gflops,"stop")
1331      write(msg,'(2(a,f8.2))')"scattered ncwrite, cpu: ",cpu,", wall: ",wall
1332      call wrtout(std_out,msg,"PERS")
1333 #endif
1334 ! HAVE_NETCDF
1335    end if !nctk_has_mpiio
1336  end if
1337 
1338  call crystal_free(crystal)
1339  if (response == 0) call ebands_free(gs_ebands)
1340 
1341  ABI_FREE(occ3d)
1342 
1343  DBG_EXIT("COLL")
1344 
1345 end subroutine cg_ncwrite

m_iowf/kg2seqblocks [ Functions ]

[ Top ] [ m_iowf ] [ Functions ]

NAME

 kg2seqblocks

FUNCTION

INPUTS

  npwtot_k
  npw_k=
  kg(3,npw_k)=reduced planewave coordinates.
  gmpi2seq
  comm_fft=FFT communicator

OUTPUT

  start_pwblock
  count_pwblock
  gblock(:,:)

PARENTS

CHILDREN

      xmpi_sum_master

SOURCE

1475 subroutine kg2seqblocks(npwtot_k,npw_k,kg_k,gmpi2seq,comm_fft,start_pwblock,count_pwblock,gblock)
1476 
1477 
1478 !This section has been created automatically by the script Abilint (TD).
1479 !Do not modify the following lines by hand.
1480 #undef ABI_FUNC
1481 #define ABI_FUNC 'kg2seqblocks'
1482 !End of the abilint section
1483 
1484  implicit none
1485 
1486 !Arguments ------------------------------------
1487 !scalars
1488  integer,intent(in) :: npwtot_k,npw_k,comm_fft
1489  integer,intent(out) :: start_pwblock,count_pwblock
1490 !arrays
1491  integer,intent(in) :: kg_k(3,npw_k),gmpi2seq(npw_k)
1492  integer,allocatable,intent(out) :: gblock(:,:)
1493 
1494 !Local variables-------------------------------
1495 !scalars
1496  integer :: me_fft,nproc_fft,rank,ig,igseq,maxnpw,ierr
1497 !arrays
1498  integer,allocatable :: igstart_rank(:),gbuf(:,:)
1499 
1500 ! *************************************************************************
1501 
1502  me_fft = xmpi_comm_rank(comm_fft); nproc_fft = xmpi_comm_size(comm_fft)
1503 
1504  ! Handle sequential case.
1505  if (nproc_fft == 1) then
1506    start_pwblock = 1; count_pwblock = npw_k
1507    ABI_MALLOC(gblock, (3, npw_k))
1508    gblock(:,:) = kg_k
1509    return ! DOH
1510  end if
1511 
1512  ABI_MALLOC(igstart_rank, (0:nproc_fft))
1513  igstart_rank = [(1 + (npwtot_k/nproc_fft)*rank, rank=0,nproc_fft-1), 1 + npwtot_k]
1514 
1515  ! Get max dimension for workspace array
1516  ! Cannot use npwtot_k / nproc_fft because G-vectors are not equally distributed.
1517  call xmpi_max(npw_k, maxnpw, comm_fft, ierr)
1518  ABI_MALLOC(gbuf, (3, maxnpw))
1519 
1520  do rank=0,nproc_fft-1
1521    start_pwblock = igstart_rank(rank)
1522    count_pwblock = igstart_rank(rank+1) - igstart_rank(rank)
1523 
1524    gbuf = 0
1525    do ig=1,npw_k
1526      igseq = gmpi2seq(ig)
1527      if (igseq >= igstart_rank(rank) .and. igseq < igstart_rank(rank+1)) then
1528        igseq = igseq - start_pwblock + 1
1529        !ABI_CHECK(igseq <= maxnpw, "boom")
1530        gbuf(:,igseq) = kg_k(:,ig)
1531      end if
1532    end do
1533    call xmpi_sum_master(gbuf,rank,comm_fft,ierr)
1534 
1535    if (me_fft == rank) then
1536      ABI_STAT_MALLOC(gblock, (3, count_pwblock), ierr)
1537      ABI_CHECK(ierr==0, "oom in gblock")
1538      gblock = gbuf(:, :count_pwblock)
1539    end if
1540  end do
1541 
1542  start_pwblock = igstart_rank(me_fft)
1543  count_pwblock = igstart_rank(me_fft+1) - igstart_rank(me_fft)
1544 
1545  ABI_FREE(gbuf)
1546  ABI_FREE(igstart_rank)
1547 
1548 end subroutine kg2seqblocks

m_iowf/ncwrite_eigen1_occ [ Functions ]

[ Top ] [ m_iowf ] [ Functions ]

NAME

 ncwrite_eigen1_occ

FUNCTION

  Write the first order DFPT eigenvalues and the occupations.

INPUTS

  ncid=Netcdf file handler.
  nband(nkpt*nsppol)=Number of bands.
  mband=maximum number of bands
  nkpt=Total number of k points
  nsppol=1 for unpolarized, 2 for spin-polarized
  eigen((2*mband**2*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point
  occ3d(mband*nkpt*nsppol)=occupations for all bands at each k point
    Note that occ3d differes from the occ arrays used in the rest of the code.
    occ3d is an array with constant stride `mband` whereas abinit (for reasons that are not clear to me)
    packs the occupations in a 1d vector with a k-dependent separation (nband_k).
    occ and occ3d differ only if nband is k-dependent but you should never assume this, hence
    remember to *convert* occ into occ3d before calling this routine.

PARENTS

      m_iowf

CHILDREN

      xmpi_sum_master

SOURCE

1379 subroutine ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d)
1380 
1381 
1382 !This section has been created automatically by the script Abilint (TD).
1383 !Do not modify the following lines by hand.
1384 #undef ABI_FUNC
1385 #define ABI_FUNC 'ncwrite_eigen1_occ'
1386 !End of the abilint section
1387 
1388  implicit none
1389 
1390 !Arguments ------------------------------------
1391 !scalars
1392  integer,intent(in) :: ncid,mband,nkpt,nsppol
1393 !arrays
1394  integer, intent(in) :: nband(nkpt*nsppol)
1395  real(dp),intent(in) :: eigen(2*mband**2*nkpt*nsppol),occ3d(mband,nkpt,nsppol)
1396 
1397 !Local variables-------------------------------
1398 !scalars
1399  integer :: idx,spin,ikpt,nband_k,ib2,ib1
1400  integer :: ncerr,occ_varid,h1mat_varid
1401 !arrays
1402  real(dp),allocatable :: h1mat(:,:,:,:,:)
1403 
1404 ! *************************************************************************
1405 
1406  ! Declare h1 array with abinit conventions
1407  ! Cannot use a 3D array since eigen1 are packed and one could have different number of bands
1408  ! per k-points. (this is not an official etsf-io variable)!
1409  ! elements in eigen are packed in the first positions.
1410  ! Remember that eigen are not MPI-distributed so no communication is needed
1411  idx=1
1412  ABI_CALLOC(h1mat, (2,mband,mband,nkpt,nsppol))
1413  do spin=1,nsppol
1414    do ikpt=1,nkpt
1415      nband_k = nband(ikpt + (spin-1)*nkpt)
1416      do ib2=1,nband_k
1417        do ib1=1,nband_k
1418           h1mat(:,ib1,ib2,ikpt,spin) = eigen(idx:idx+1)
1419           idx=idx+2
1420         end do
1421      end do
1422    end do
1423  end do
1424 
1425  !write(std_out,*)"in occ1"
1426  NCF_CHECK(nctk_set_defmode(ncid))
1427 
1428  ncerr = nctk_def_arrays(ncid, nctkarr_t('h1_matrix_elements', "dp", &
1429 &"complex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins"))
1430  NCF_CHECK(ncerr)
1431 
1432  NCF_CHECK(nf90_inq_varid(ncid, "occupations", occ_varid))
1433  NCF_CHECK(nf90_inq_varid(ncid, "h1_matrix_elements", h1mat_varid))
1434 
1435  ! Write data
1436  NCF_CHECK(nctk_set_datamode(ncid))
1437  NCF_CHECK_MSG(nf90_put_var(ncid, occ_varid, occ3d), "putting occ3d")
1438  NCF_CHECK_MSG(nf90_put_var(ncid, h1mat_varid, h1mat), "putting h1mat")
1439 
1440  ABI_FREE(h1mat)
1441 
1442 end subroutine ncwrite_eigen1_occ

m_iowf/outwf [ Functions ]

[ Top ] [ m_iowf ] [ Functions ]

NAME

 outwf

FUNCTION

 Conduct output of a "wave-functions" file.
  - Compute the maximal residual
  - Then open a permanent file wff2 for final output of wf data
  - Create a new header for the file.
  - Write wave-functions (and energies)

INPUTS

  cg(2,mcg)=wavefunction array (storage if nkpt>1)
  dtset <type(dataset_type)>=all input variables for this dataset
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  eigen( (2*mband)**response *mband*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point
  filnam= character string giving the root to form the name of the
   output WFK or WFQ file if response==0, otherwise it is the filename.
  hdr <type(hdr_type)>=the header of wf, den and pot files
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kptns(3,nkpt)=k points in terms of recip primitive translations
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mkmem=Number of k-points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpw=maximum number of plane waves
  natom=number of atoms in unit cell
  nband=number of bands
  nkpt=number of k points
  npwarr(nkpt)=number of planewaves in basis and on boundary for each k
  nsppol=1 for unpolarized, 2 for spin-polarized
  nstep=desired number of electron iteration steps
  occ(mband*nkpt*nsppol)=occupations for all bands at each k point
  resid(mband*nkpt*nsppol)=squared residuals for each band and k point
   where resid(n,k)=|<C(n,k)|(H-e(n,k))|C(n,k)>|^2 for the ground state
  response: if == 0, GS wavefunctions , if == 1, RF wavefunctions
  unwff2=unit for output of wavefunction
  wfs <type(wvl_projector_type)>=wavefunctions information for wavelets.

OUTPUT

  (only writing)

NOTES

 * The name of the file wff2 might be the same as that of the file wff1.

PARENTS

      berryphase_new,dfpt_looppert,gstate

CHILDREN

      xmpi_sum_master

SOURCE

118 subroutine outwf(cg,dtset,psps,eigen,filnam,hdr,kg,kptns,mband,mcg,mkmem,&
119 &                mpi_enreg,mpw,natom,nband,nkpt,npwarr,&
120 &                nsppol,occ,resid,response,unwff2,&
121 &                wfs,wvl)
122 
123 
124 !This section has been created automatically by the script Abilint (TD).
125 !Do not modify the following lines by hand.
126 #undef ABI_FUNC
127 #define ABI_FUNC 'outwf'
128 !End of the abilint section
129 
130  implicit none
131 
132 !Arguments ------------------------------------
133 !scalars
134  integer,intent(in) :: mband,mcg,mkmem,mpw,natom,nkpt,nsppol,response,unwff2
135 !integer,intent(in) :: nstep
136  character(len=*),intent(in) :: filnam
137  type(MPI_type),intent(in) :: mpi_enreg
138  type(dataset_type),intent(in) :: dtset
139  type(pseudopotential_type),intent(in) :: psps
140  type(hdr_type), intent(inout) :: hdr
141  type(wvl_wf_type),intent(in) :: wfs
142  type(wvl_internal_type), intent(in) :: wvl
143 !arrays
144  integer, intent(in) :: kg(3,mpw*mkmem),nband(nkpt*nsppol),npwarr(nkpt)
145  real(dp), intent(inout) :: cg(2,mcg)
146  real(dp), intent(in) :: eigen((2*mband)**response*mband*nkpt*nsppol),kptns(3,nkpt)
147  real(dp), intent(in) :: occ(mband*nkpt*nsppol),resid(mband*nkpt*nsppol)
148 
149 !Local variables-------------------------------
150  integer,parameter :: nkpt_max=50
151  integer :: iomode,action,band_index,fform,formeig,iband,ibdkpt,icg,iat,iproj
152  integer :: ierr,ii,ikg,ikpt,spin,master,mcg_disk,me,me0,my_nspinor
153  integer :: nband_k,nkpt_eff,nmaster,npw_k,option,rdwr,sender,source !npwtot_k,
154  integer :: spaceComm,spaceComm_io,spacecomsender,spaceWorld,sread,sskip,tim_rwwf,xfdim2
155 #ifdef HAVE_MPI
156  integer :: ipwnbd
157 #endif
158  real(dp) :: residk,residm,resims,cpu,wall,gflops
159  logical :: ihave_data,iwrite,iam_master,done
160  character(len=500) :: msg
161  type(wffile_type) :: wff2
162  character(len=fnlen) :: path
163 !arrays
164  integer,allocatable :: kg_disk(:,:)
165  real(dp) :: tsec(2)
166  real(dp),allocatable :: cg_disk(:,:),eig_k(:),occ_k(:)
167 #ifdef HAVE_NETCDF
168  integer :: ncid, ncerr, kg_varid, mpw_disk, npwk_disk, timrev
169  real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:),vkbsign(:,:)
170  type(crystal_t) :: crystal
171 #endif
172 
173 ! *************************************************************************
174 !For readability of the source file, define a "me" variable also in the sequential case
175 
176  DBG_ENTER("COLL")
177 
178  xfdim2 = natom+4
179 !Init mpi_comm
180  spaceWorld= mpi_enreg%comm_cell
181  spaceComm=spaceWorld
182  spaceComm_io=xmpi_comm_self
183 
184  if (mpi_enreg%paral_kgb==1 ) spaceComm_io= mpi_enreg%comm_bandspinorfft
185  if (mpi_enreg%paral_kgb==1 ) spaceComm= mpi_enreg%comm_cell
186 
187 !Paral_kgb=1 and Fortran-I/O is not supported (only for testing purpose)
188  if (mpi_enreg%paral_kgb==1.and.dtset%iomode==IO_MODE_FORTRAN) then
189    spaceWorld=mpi_enreg%comm_kpt
190    write(msg,'(7a)') &
191 &   'WF file is written using standard Fortran I/O',ch10,&
192 &   'and Kpt-band-FFT parallelization is active !',ch10,&
193 &   'This is only allowed for testing purposes.',ch10,&
194 &   'The produced WF file will be incomplete and not useable.'
195    MSG_WARNING(msg)
196  end if
197 
198 !If parallel HF calculation
199  if (mpi_enreg%paral_hf==1 ) spaceComm_io= mpi_enreg%comm_hf
200  if (mpi_enreg%paral_hf==1 ) spaceComm= mpi_enreg%comm_cell
201 
202 !Paral_hf=1 and Fortran-I/O is not supported (copy from paral_kgb... not tested)
203  if (mpi_enreg%paral_hf==1.and.dtset%iomode==IO_MODE_FORTRAN) then
204    spaceWorld=mpi_enreg%comm_kpt
205    write(msg,'(7a)') &
206 &   'WF file is written using standard Fortran I/O',ch10,&
207 &   'and HF parallelization is active !',ch10,&
208 &   'This is only allowed for testing purposes.',ch10,&
209 &   'The produced WF file will be incomplete and not useable.'
210    MSG_WARNING(msg)
211  end if
212 
213  ! check consistency between dimensions and input hdr.
214  !ABI_CHECK(mband == maxval(hdr%nband), "hdr:mband")
215  !ABI_CHECK(nkpt == hdr%nkpt, "hdr:nkpt")
216  !ABI_CHECK(nsppol == hdr%nsppol, "hdr:nsppol")
217  !ABI_CHECK(all(hdr%npwarr == npwarr), "hdr:npwarr")
218  !ABI_CHECK(all(hdr%nband == nband), "hdr:nband")
219  !ABI_CHECK(maxval(hdr%npwarr) == mpw, "hdr:nband")
220 
221 !Init me
222  me=mpi_enreg%me_kpt
223  me0=me
224 !Define master
225  master=0
226 
227  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
228  tim_rwwf =0
229  source = master
230  sread = master
231  iam_master=(master==me)
232  iwrite=iam_master
233  sender=-1
234 
235 !Compute mean square and maximum residual over all bands and k points and spins
236 !(disregard k point weights and occupation numbers here)
237  band_index=sum(nband(1:nkpt*nsppol))
238  resims=sum(resid(1:band_index))/dble(band_index)
239 
240 !Find largest residual over bands, k points, and spins, except for nbdbuf highest bands
241 !Already AVAILABLE in hdr ?!
242  ibdkpt=1
243  residm=zero
244  do spin=1,nsppol
245    do ikpt=1,nkpt
246      nband_k=nband(ikpt+(spin-1)*nkpt)
247      nband_k=max(1,nband_k-dtset%nbdbuf)
248      residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_k-1)))
249      ibdkpt=ibdkpt+nband_k
250    end do
251  end do
252 
253  write(msg,'(a,2p,e12.4,a,e12.4)')' Mean square residual over all n,k,spin= ',resims,'; max=',residm
254  call wrtout(ab_out,msg,'COLL')
255  call wrtout(std_out,msg,'COLL')
256 
257  band_index=0
258  nkpt_eff=nkpt
259  if( (dtset%prtvol==0 .or. dtset%prtvol==1) .and. nkpt_eff>nkpt_max ) nkpt_eff=nkpt_max
260 
261 !Loop over spin again
262  do spin=1,nsppol
263 !  Give (squared) residuals for all bands at each k
264    do ikpt=1,nkpt
265      nband_k=nband(ikpt+(spin-1)*nkpt)
266 !    Will not print all residuals when prtvol=0 or 1
267      if(ikpt<=nkpt_eff)then
268 !      Find largest residual over all bands for given k point
269        residk=maxval(resid(1+band_index:nband_k+band_index))
270        write(msg,'(1x,3f8.4,3x,i2,1p,e13.5,a)')kptns(1:3,ikpt),spin,residk,' kpt; spin; max resid(k); each band:'
271        if(dtset%prtvol>=2)then
272          call wrtout(ab_out,msg,'COLL')
273        endif
274        call wrtout(std_out,msg,'COLL')
275        do ii=0,(nband_k-1)/8
276          write(msg,'(1x,1p,8e9.2)')(resid(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8))
277          if(dtset%prtvol>=2)then
278            call wrtout(ab_out,msg,'COLL')
279          endif
280          call wrtout(std_out,msg,'COLL')
281        end do
282      else if(ikpt==nkpt_eff+1)then
283        write(msg,'(2a)')' outwf : prtvol=0 or 1, do not print more k-points.',ch10
284        if(dtset%prtvol>=2)then
285          call wrtout(ab_out,msg,'COLL')
286        endif
287        call wrtout(std_out,msg,'COLL')
288      end if
289      band_index=band_index+nband_k
290    end do
291  end do
292 
293 !Will write the wavefunction file only when nstep>0
294 !MT 07 2015: writing reactivated when nstep=0
295 !if (nstep>0 .and. dtset%prtwf/=0) then
296  if (dtset%prtwf/=0) then
297 
298    ! Only the master write the file, except if MPI I/O, but the
299    ! full wff dataset should be provided to WffOpen in this case
300    iomode=IO_MODE_FORTRAN_MASTER
301    if (dtset%iomode==IO_MODE_MPI)  iomode = IO_MODE_MPI
302    if (dtset%iomode==IO_MODE_ETSF) iomode = IO_MODE_ETSF
303 
304 #ifdef HAVE_NETCDF
305    if (dtset%iomode == IO_MODE_ETSF .and. dtset%usewvl == 0) then
306      call cg_ncwrite(filnam,hdr,dtset,response,mpw,mband,nband,nkpt,nsppol,&
307        dtset%nspinor,mcg,mkmem,eigen,occ,cg,npwarr,kg,mpi_enreg,done)
308 #ifdef HAVE_NETCDF_DEFAULT
309      ABI_CHECK(done, "cg_ncwrite must handle the output of the WFK file.")
310 #endif
311 
312      ! Write KB form factors. Only master works. G-vectors are read from file to avoid
313      ! having to deal with paral_kgb distribution.
314      if (me == master .and. dtset%prtkbff == 1 .and. dtset%iomode == IO_MODE_ETSF .and. dtset%usepaw == 0) then
315        ABI_CHECK(done, "cg_ncwrite was not able to generate WFK.nc in parallel. Perphase hdf5 is not working")
316        path = nctk_ncify(filnam)
317        call wrtout(std_out, sjoin("Writing KB form factors to:", path))
318        NCF_CHECK(nctk_open_modify(ncid, path, xmpi_comm_self))
319        NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid))
320        mpw_disk = maxval(hdr%npwarr)
321 
322        ! Dimensions needed by client code to allocate memory when reading.
323        ncerr = nctk_def_dims(ncid, [ &
324          nctkdim_t("mproj", psps%mproj), &
325          nctkdim_t("mpsang", psps%mpsang), &
326          nctkdim_t("mpssoang", psps%mpssoang), &
327          nctkdim_t("lnmax", psps%lnmax), &
328          nctkdim_t("lmnmax", psps%lmnmax) &
329        ])
330        NCF_CHECK(ncerr)
331 
332        ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "mpspso"])
333        NCF_CHECK(ncerr)
334 
335        ! Write indlmn table (needed to access vkb arrays)
336        ncerr = nctk_def_arrays(ncid, [ &
337          nctkarr_t("indlmn", "int", "six, lmnmax, number_of_atom_species"), &
338          nctkarr_t("vkbsign", "dp", "lnmax, number_of_atom_species"), &
339          nctkarr_t("vkb", "dp", "max_number_of_coefficients, lnmax, number_of_atom_species, number_of_kpoints"), &
340          nctkarr_t("vkbd", "dp", "max_number_of_coefficients, lnmax, number_of_atom_species, number_of_kpoints") &
341        ], defmode=.True.)
342        NCF_CHECK(ncerr)
343 
344        ! Switch to write mode.
345        NCF_CHECK(nctk_set_datamode(ncid))
346        NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "indlmn"), psps%indlmn))
347 
348        ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
349          "mpspso"], &
350          [psps%mpspso])
351        NCF_CHECK(ncerr)
352 
353        ! Calculate KB form factors and derivatives.
354        ! The arrays are allocated with lnmax to support pseudos with more than projector.
355        ! Note that lnmax takes into account lloc hence arrays are in packed form and one should be
356        ! accessed with the indices provided by psps%indlmn.
357        ABI_MALLOC(vkbsign, (psps%lnmax, psps%ntypat))
358        ABI_MALLOC(vkb, (mpw_disk, psps%lnmax, psps%ntypat))
359        ABI_MALLOC(vkbd, (mpw_disk, psps%lnmax, psps%ntypat))
360        ABI_MALLOC(kg_disk, (3, mpw_disk))
361 
362        timrev = 2 ! FIXME: Use abinit convention for timrev
363        call crystal_from_hdr(crystal, hdr, timrev)
364 
365        ! For each k-point: read full G-vector list from file, compute KB data and write to file.
366        do ikpt=1,nkpt
367          npwk_disk = hdr%npwarr(ikpt)
368          NCF_CHECK(nf90_get_var(ncid, kg_varid, kg_disk, start=[1, 1, ikpt], count=[3, npwk_disk, 1]))
369          vkb = zero; vkbd = zero
370          call calc_vkb(crystal, psps, kptns(:, ikpt), npwk_disk, mpw_disk, kg_disk, vkbsign, vkb, vkbd)
371 
372          if (ikpt == 1) then
373            ! This for the automatic tests.
374            NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkbsign"), vkbsign))
375            if(dtset%prtvol>=2)then
376              write(msg,'(a)') 'prtkbff: writing first and last G-components of the KB form factors'
377              call wrtout(ab_out,msg,'COLL')
378              do iat=1,psps%ntypat
379                write(msg,'(a10,i5)') 'atype ',iat
380                call wrtout(ab_out,msg,'COLL')
381                do iproj=1,psps%lnmax
382                  write(msg,'(a10,i5,a,a10,e12.4,a,3(a10,2e12.4,a))') &
383                         'projector ', iproj,ch10, &
384                         'vkbsign   ', vkbsign(iproj,iat), ch10, &
385                         'vkb       ', vkb(1,iproj,iat),  vkb(npwk_disk,:,iat), ch10, &
386                         'vkbd      ', vkbd(1,iproj,iat), vkbd(npwk_disk,:,iat), ''
387                  call wrtout(ab_out,msg,'COLL')
388                end do
389              end do
390            end if
391          end if
392 
393          NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkb"), vkb, start=[1, 1, 1, ikpt]))
394          NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkbd"), vkbd, start=[1, 1, 1, ikpt]))
395        end do
396        NCF_CHECK(nf90_close(ncid))
397 
398        ABI_FREE(kg_disk)
399        ABI_FREE(vkbsign)
400        ABI_FREE(vkb)
401        ABI_FREE(vkbd)
402        call crystal_free(crystal)
403      end if
404 
405      if (done) return
406      ! If cg_ncwrite cannot handle the IO because HDF5 + MPI-IO support is missing, we fallback to Fortran + MPI-IO.
407      msg = "Could not produce a netcdf file in parallel (MPI-IO support is missing). Will fallback to MPI-IO with Fortran"
408      MSG_WARNING(msg)
409      iomode=IO_MODE_MPI
410    end if
411 #endif
412 
413    call cwtime(cpu, wall, gflops, "start")
414 
415    write(msg,'(4a,i0)')ch10,' outwf: write wavefunction to file ',trim(filnam),", with iomode ",iomode
416    call wrtout(std_out,msg,'COLL')
417 
418    ! Create an ETSF file for the wavefunctions
419    if (iomode == IO_MODE_ETSF) then
420      ABI_CHECK(xmpi_comm_size(spaceComm) == 1, "Legacy etsf-io code does not support nprocs > 1")
421 #ifdef HAVE_ETSF_IO
422      call abi_etsf_init(dtset, filnam, 2, .true., hdr%lmn_size, psps, wfs)
423      !call crystal_from_hdr(crystal, hdr, 2)
424      !NCF_CHECK(crystal_ncwrite_path(crystal, nctk_ncify(filnam)))
425      !call crystal_free(crystal)
426      !ncerr = ebands_ncwrite_path(gs_ebands, filname, ncid)
427      !NCF_CHECK(ncerr)
428 #else
429      MSG_ERROR("ETSF_IO is not activated")
430      ABI_UNUSED(psps%ntypat)
431 #endif
432    end if
433 
434    call WffOpen(iomode,spaceComm,filnam,ierr,wff2,master,me0,unwff2,spaceComm_io)
435 !  Conduct wavefunction output to wff2
436 
437    ABI_ALLOCATE(kg_disk,(3,mpw))
438 
439    mcg_disk=mpw*my_nspinor*mband
440    formeig=0; if (response==1) formeig=1
441 
442    ABI_ALLOCATE(eig_k,( (2*mband)**formeig * mband))
443    ABI_ALLOCATE(occ_k,(mband))
444 
445 #ifdef HAVE_MPI
446    call xmpi_barrier(spaceComm)
447 !  Compute mband and mpw
448    ABI_STAT_ALLOCATE(cg_disk,(2,mcg_disk), ierr)
449    ABI_CHECK(ierr==0, "out of memory in cg_disk")
450 #endif
451 
452    band_index=0
453    icg=0
454    if(mpi_enreg%paralbd==0) tim_rwwf=6
455    if(mpi_enreg%paralbd==1) tim_rwwf=12
456 
457 !  Write header info for new wf file
458    rdwr=2
459    if (dtset%usewvl==0) then
460      fform=2
461    else
462      fform = 200 ! Use 200 as radical for naming file format used by wavelets.
463    end if
464 
465    if (wff2%iomode < 2) then
466      call hdr_io(fform,hdr,rdwr,wff2)
467      call WffKg(wff2,1)
468    else if (wff2%iomode==IO_MODE_ETSF .and. iam_master) then
469 #ifdef HAVE_NETCDF
470      NCF_CHECK(hdr_ncwrite(hdr, wff2%unwff, fform, nc_define=.True.))
471 #endif
472    end if
473 
474    do spin=1,nsppol
475      ikg=0
476 
477      do ikpt=1,nkpt
478        nband_k=nband(ikpt+(spin-1)*nkpt)
479        npw_k=npwarr(ikpt)
480 
481 #ifdef HAVE_MPI
482        if (dtset%usewvl == 0) then
483          call xmpi_barrier(spaceWorld)
484 
485 !        Must transfer the wavefunctions to the master processor
486 !        Separate sections for paralbd=1 or other values ; might be merged
487          if(mpi_enreg%paralbd==0)then
488            nmaster=0
489            source=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,spin))
490            ihave_data=.false.
491            if(source==me)ihave_data=.true.
492            action=0
493 !          I am the master node, and I have the data in cg or cg_disk
494            if((iam_master).and.(ihave_data))action=1
495 !          I am not the master, and I have the data => send to master
496            if((.not.iam_master).and.(ihave_data))action=2
497 !          I am the master, and I receive the data
498            if((iam_master).and.(.not.ihave_data))action=3
499 
500 !          I have the data in cg or cg_disk ( MPI_IO case)
501            if (iomode==IO_MODE_MPI) then
502              action = 0
503              sender=-1
504              iwrite=.false.
505              if (ihave_data)then
506                action=1
507                iwrite=.true.
508                sender=me
509              end if
510            end if
511 
512 !          I am the master node, and I have the data in cg or cg_disk
513 !          I have the data in cg or cg_disk ( MPI_IO case)
514            if(action==1)then
515 !            Copy from kg to kg_disk
516              kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
517 !            Copy from cg to cg_disk
518              do ipwnbd=1,nband_k*npw_k*my_nspinor
519                cg_disk(1,ipwnbd)=cg(1,ipwnbd+icg)
520                cg_disk(2,ipwnbd)=cg(2,ipwnbd+icg)
521              end do
522            end if
523 
524 !          I am not the master, and I have the data => send to master
525 !          I am the master, and I receive the data
526            if ( action==2.or.action==3) then
527              !write(std_out,*)npw_k,nband_k
528              call timab(48,1,tsec)
529              if(action==2)then
530                call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,source,kg_disk,nmaster,spaceWorld,ierr)
531                call xmpi_exch(cg(:,icg+1:icg+nband_k*npw_k*my_nspinor),2*nband_k*npw_k*my_nspinor, &
532 &               source,cg_disk,nmaster,spaceWorld,ierr)
533              else
534                call xmpi_exch(kg_disk,3*npw_k,source,kg_disk,nmaster,spaceWorld,ierr)
535                call xmpi_exch(cg_disk,2*nband_k*npw_k*my_nspinor,source,cg_disk,nmaster,spaceWorld,ierr)
536              end if
537              call timab(48,2,tsec)
538            end if
539 
540 
541          else if(mpi_enreg%paralbd==1)then
542            nmaster=0
543 #ifdef HAVE_MPI_IO
544            sender=IO_MODE_FORTRAN_MASTER
545            if( iomode==IO_MODE_MPI) then
546              nmaster=mpi_enreg%proc_distrb(ikpt,1,spin)
547              sender=nmaster
548            end if
549 #endif
550 
551 !          Note the loop over bands
552            do iband=1,nband_k
553 
554 !            The message passing related to kg is counted as one band
555              action=0
556 
557 !            I am the master node, and I have the data in cg or cg_disk
558              if( mpi_enreg%proc_distrb(ikpt,iband,spin)==nmaster .and. me==nmaster) then
559                action=1
560 !              I am not the master, and I have the data => send to master
561              elseif( mpi_enreg%proc_distrb(ikpt,iband,spin)==me .and. me/=nmaster ) then
562                action = 2
563 !              I am the master, and I receive the data
564              elseif( mpi_enreg%proc_distrb(ikpt,iband,spin)/=me .and. me==nmaster ) then
565                action=3
566              end if
567 
568              if(action==1) then
569 !              I am the master node, and I have the data in cg or cg_disk
570 !              Copy from kg to kg_disk
571                if(iband==1)kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
572 !              Copy from cg to cg_disk
573                do ipwnbd=1,npw_k*my_nspinor
574                  cg_disk(1,(iband-1)*npw_k*my_nspinor+ipwnbd) = cg(1,(iband-1)*npw_k*my_nspinor+ipwnbd+icg)
575                  cg_disk(2,(iband-1)*npw_k*my_nspinor+ipwnbd) = cg(2,(iband-1)*npw_k*my_nspinor+ipwnbd+icg)
576                end do
577              end if  ! action=1
578 
579              if ( action==2.or.action==3) then
580 !              action=2 :  I am not the master, and I have the data => send to master
581 !              action=3 :  I am the master, and I receive the data
582                call timab(48,1,tsec)
583                if ( iband == 1 ) then
584                  if (action==2) then
585                    call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,mpi_enreg%proc_distrb(ikpt,iband,spin), &
586 &                   kg_disk,nmaster,spaceWorld,ierr)
587                  else
588                    call xmpi_exch(kg_disk,3*npw_k,mpi_enreg%proc_distrb(ikpt,iband,spin),  &
589 &                   kg_disk,nmaster,spaceWorld,ierr)
590                  end if
591                end if       ! iband =1
592                ipwnbd=(iband-1)*npw_k*my_nspinor
593                if (action==2) then
594                  call xmpi_exch( cg(:,ipwnbd+icg+1:ipwnbd+icg+npw_k*my_nspinor),2*npw_k*my_nspinor &
595 &                 ,mpi_enreg%proc_distrb(ikpt,iband,spin)                    &
596 &                 ,cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),nmaster,spaceWorld,ierr)
597                else
598                  call xmpi_exch( cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),2*npw_k*my_nspinor    &
599 &                 ,mpi_enreg%proc_distrb(ikpt,iband,spin)                    &
600 &                 ,cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),nmaster,spaceWorld,ierr)
601                end if
602 
603                call timab(48,2,tsec)
604              end if        ! action=2 or action=3
605 
606              if(iomode==IO_MODE_MPI) then
607 !              I have the data in cg or cg_disk
608                iwrite=.false.
609                if (nmaster == me) iwrite=.true.
610              end if
611 
612            end do ! End of loop over bands
613          end if ! End of paralbd=1
614        end if
615 #endif
616 
617 !      Only the master will write to disk the final output wf file.
618 !      in MPI_IO case only iwrite will write to disk the final output wf file.
619        if(iwrite) then
620 !        write(std_out,*) 'outwf : I am master and will write wf file'
621          if(formeig==0)then
622            eig_k(1:nband_k)=eigen(1+band_index:nband_k+band_index)
623            occ_k(1:nband_k)=occ(1+band_index:nband_k+band_index)
624          else
625            eig_k(1:2*nband_k*nband_k)=eigen(1+band_index:2*nband_k*nband_k+band_index)
626          end if
627          option=2
628          if(dtset%prtwf==3)option=5
629 !        if (dtset%prtwf == 2 .and. mkmem/=0) option=4
630 
631          if (dtset%usewvl == 0) then
632 #ifdef HAVE_MPI
633            call rwwf(cg_disk,eig_k,formeig,0,0,ikpt,spin,kg_disk,mband,mcg_disk,mpi_enreg, &
634 &           nband_k, nband_k,npw_k,my_nspinor,occ_k,option,1,tim_rwwf,wff2)
635 
636 #else
637            kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
638            call rwwf(cg,eig_k,formeig,0,icg,ikpt,spin,kg_disk,mband,mcg,mpi_enreg,nband_k, &
639 &           nband_k, npw_k,my_nspinor,occ_k,option,1,tim_rwwf,wff2)
640 #endif
641          else
642            call wvl_write(dtset,eigen,mpi_enreg,option,hdr%rprimd,wff2,wfs,wvl,hdr%xred)
643          end if
644        end if
645 
646 !      The wavefunctions for the present k point and spin are written
647        if(response==0)band_index=band_index+nband_k
648        if(response==1)band_index=band_index+2*nband_k*nband_k
649 
650        sskip=1
651 #ifdef HAVE_MPI
652        if (dtset%usewvl == 0) then
653          sskip=0
654          if(.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me)))sskip=1
655        end if
656 #endif
657        if(sskip==1)then
658          icg=icg+npw_k*my_nspinor*nband_k
659          ikg=ikg+npw_k
660        end if
661 
662 
663 #ifdef HAVE_MPI_IO
664        spacecomsender=spaceComm
665        if (mpi_enreg%paral_kgb==1) spacecomsender =mpi_enreg%comm_kpt
666        if (mpi_enreg%paral_hf==1) spacecomsender =mpi_enreg%comm_kpt
667        call WffOffset(wff2,sender,spacecomsender,ierr)
668 #endif
669 
670      end do ! ikpt
671    end do ! spin
672 
673    ABI_DEALLOCATE(kg_disk)
674 #ifdef HAVE_MPI
675    ABI_DEALLOCATE(cg_disk)
676 #endif
677 
678    ABI_DEALLOCATE(eig_k)
679    ABI_DEALLOCATE(occ_k)
680 
681 !  Close the wavefunction file (and do NOT delete it !)
682    !if (wff2%iomode /= IO_MODE_NETCDF) then
683    call WffClose(wff2,ierr)
684    !end if
685 
686    call cwtime(cpu,wall,gflops,"stop")
687    write(msg,'(a,i0,2(a,f8.2),a)')" outwf with iomode: ",iomode,", cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
688    call wrtout(std_out,msg,"PERS")
689  end if ! End condition of nstep>0
690 
691  ! Block here because we might need to read the WFK file in the caller.
692  call xmpi_barrier(mpi_enreg%comm_cell)
693 
694  DBG_EXIT("COLL")
695 
696 end subroutine outwf