TABLE OF CONTENTS


ABINIT/m_ioarr [ Modules ]

[ Top ] [ Modules ]

NAME

 m_ioarr

FUNCTION

  This module provides routines to read/write arrays given on the FFT mesh (densities, potentials ...).
  The code supports both Fortran files as well as netcdf files in a transparent way.
  The appropriate IO layer is selected from files extensions: netcdf primitives are used if the
  file ends with `.nc`. If all the other cases we read/write files in Fortran format.
  MPI-IO primitives are used when the FFT arrays are MPI distributed.

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, MVer, MT, 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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 MODULE m_ioarr
27 
28  use defs_basis
29  use m_abicore
30  use m_xmpi
31  use m_wffile
32  use m_errors
33  use m_nctk
34  use m_dtset
35  use m_crystal
36  use m_ebands
37  use m_hdr
38  use m_pawrhoij
39 #ifdef HAVE_MPI2
40  use mpi
41 #endif
42  use netcdf
43 
44  use defs_abitypes,   only : mpi_type
45  use defs_datatypes,  only : ebands_t
46  use defs_wvltypes,   only : wvl_denspot_type
47  use m_time,          only : cwtime, cwtime_report, timab
48  use m_io_tools,      only : iomode_from_fname, iomode2str, open_file, get_unit
49  use m_fstrings,      only : sjoin, itoa, endswith, ltoa
50  use m_numeric_tools, only : interpolate_denpot
51  use m_geometry,      only : metric
52  use m_mpinfo,        only : destroy_mpi_enreg, ptabs_fourdp, initmpi_seq
53  use m_distribfft,    only : init_distribfft_seq
54  use m_fourier_interpol,only : fourier_interpol
55 
56  implicit none
57 
58 #ifdef HAVE_MPI1
59  include 'mpif.h'
60 #endif
61 
62  private
63 
64  public :: ioarr                     ! Read or write rho(r) or v(r), either ground-state or response-functions.
65  public :: fftdatar_write            ! Write an array in real space. IO library is automatically selected
66                                      ! from the file extension and the number of FFT processors:
67  public :: fftdatar_write_from_hdr   ! Write an array in real-space to file plus crystal_t and ebands_t
68  public :: read_rhor                 ! Read rhor from DEN file.
69  public :: fort_denpot_skip          ! Skip the header and the DEN/POT records (Fortran format)
70 
71  private :: denpot_spin_convert      ! Convert a density/potential from a spin representation to another
72 
73 CONTAINS  !====================================================================================================

m_ioarr/denpot_spin_convert [ Functions ]

[ Top ] [ m_ioarr ] [ Functions ]

NAME

  denpot_spin_convert

FUNCTION

  Convert a density/potential from a spin representation to another

INPUTS

  denpot_in(:,nspden_in)=input density//potential
  nspden_in=number of spin-component of the input density/potential
  fform=file format (density or potential)
  [istart_in]= --optional-- starting index in the denpot_in array; default is 1
  [istart_out]= --optional-- starting index in the denpot_out array; default is 1
  [nelem]= --optional-- number of elements to copy from denpot_in to denpot_out; default is all

OUTPUT

  denpot_out(:,nspden_out)=output density//potential
  nspden_out=number of spin-component of the output density/potential

NOTES

  More explicitely:
    We copy denpot_in(istar_in+1:istart_in+nelem,:)
       into denpot_out(istart_out+1:istart_out+nelem,:)

SOURCE

1311 subroutine denpot_spin_convert(denpot_in,nspden_in,denpot_out,nspden_out,fform,&
1312 &                              istart_in,istart_out,nelem) ! optional arguments
1313 
1314 !Arguments ------------------------------------
1315 !scalars
1316  integer,intent(in) :: nspden_in,nspden_out,fform
1317  integer,intent(in),optional :: istart_in,istart_out,nelem
1318 !arrays
1319  real(dp),intent(in) :: denpot_in(:,:)
1320  real(dp),intent(out) :: denpot_out(:,:)
1321 
1322 !Local variables-------------------------------
1323  integer :: iend_in,iend_out,ispden,my_istart_in,my_istart_out,my_nelem
1324  character(len=500) :: msg
1325 
1326 ! *********************************************************************
1327 
1328 !Optional arguments
1329  my_istart_in=1;if (present(istart_in)) my_istart_in=istart_in
1330  my_istart_out=1;if (present(istart_out)) my_istart_out=istart_out
1331  iend_in=size(denpot_in,1) ; iend_out=size(denpot_out,1)
1332  my_nelem=min(iend_in-my_istart_in+1,iend_out-my_istart_out+1)
1333  if (present(nelem)) my_nelem=nelem
1334 
1335 !Checks
1336  if (size(denpot_in,2)/=nspden_in) then
1337    msg='size(denpot_in,2)/=nspden_in!'
1338    ABI_BUG(msg)
1339  end if
1340  if (size(denpot_out,2)/=nspden_out) then
1341    msg='size(denpot_out,2)/=nspden_out!'
1342    ABI_BUG(msg)
1343  end if
1344  if (my_istart_in+my_nelem-1>size(denpot_in,1)) then
1345    msg='istart_in+nelem>size(denpot_in,1)!'
1346    ABI_BUG(msg)
1347  end if
1348  if (my_istart_out+my_nelem-1>size(denpot_out,1)) then
1349    msg='istart_out+nelem>size(denpot_out,1)!'
1350    ABI_BUG(msg)
1351  end if
1352 
1353 !Simple copy if the number of spin-components is unchanged...
1354  if (nspden_in==nspden_out) then
1355    do ispden=1,nspden_in
1356      denpot_out(my_istart_out:my_istart_out+my_nelem-1,ispden)= &
1357 &      denpot_in(my_istart_in:my_istart_in+my_nelem-1,ispden)
1358    end do
1359    return
1360  end if
1361 
1362 !...otherwise, we need to convert.
1363  if ((fform-1)/2==25) then
1364 
1365 !  First case: DENSITY
1366 
1367    if      (nspden_in==1.and.nspden_out==2) then
1368      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1369      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)*half
1370    else if (nspden_in==1.and.nspden_out==4) then
1371      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1372      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=zero
1373      denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1374      denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1375    else if (nspden_in==2.and.nspden_out==1) then
1376      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1377    else if (nspden_in==2.and.nspden_out==4) then
1378      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1379      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=zero
1380      denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1381      denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*two &
1382 &                                                        -denpot_in(my_istart_in:my_istart_in+my_nelem,1)
1383    else if (nspden_in==4.and.nspden_out==1) then
1384      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1385    else if (nspden_in==4.and.nspden_out==2) then
1386      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1387      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)*half &
1388 &                                                        +denpot_in(my_istart_in:my_istart_in+my_nelem-1,4)*half
1389    end if
1390 
1391  else
1392 
1393 !  Second case: POTENTIAL
1394 
1395    if      (nspden_in==1.and.nspden_out==2) then
1396      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1397      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1398    else if (nspden_in==1.and.nspden_out==4) then
1399      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1400      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1401      denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1402      denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1403    else if (nspden_in==2.and.nspden_out==1) then
1404      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)*half &
1405 &                                                        +denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*half
1406    else if (nspden_in==2.and.nspden_out==4) then
1407      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1408      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)
1409      denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1410      denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1411    else if (nspden_in==4.and.nspden_out==1) then
1412      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)*half &
1413 &                                                        +denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*half
1414    else if (nspden_in==4.and.nspden_out==2) then
1415      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1416      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)
1417    end if
1418 
1419  end if
1420 
1421 end subroutine denpot_spin_convert

m_ioarr/fftdatar_write [ Functions ]

[ Top ] [ m_ioarr ] [ Functions ]

NAME

 fftdatar_write

FUNCTION

 Write an array in real space on the FFT box to file.
 The array can be real or complex depending on cplex
 IO library is automatically selected from the file extension and the number of FFT processors:

   1) If path ends with ".nc", the netcdf library is used else Fortran format.

   2) If nproc_fft > 1, parallel IO is used (if available)

INPUTS

 varname=Name of the variable to write (used if ETSF-IO).
 path=File name
 iomode=
 hdr <type(hdr_type)>=the header of wf, den and pot files
 crystal<crystal_t>= data type gathering info on symmetries and unit cell (used if etsf_io)
 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 cplex=1 for real array, 2 for complex
 nfft=Number of FFT points treated by this node.
 nspden=Number of spin-density components.
 datar(cplex*nfft,nspden)=array on the real space FFT grid.
 mpi_enreg=information about MPI parallelization
 [ebands]<ebands_t>=data type with energies and occupations (used if etsf_io)

OUTPUT

  Only writing

NOTES

   The string passed to fftdatar_write (first argument) gives the name used to store the data in the netcdf file
   The function varname_from_fname defined in the module m_hdr.F90 gives the mapping between the Abinit
   file extension and the netcdf name e.g. foo_VHXC.nc --> vxc
   This function is used in cut3d so that we can immediately select the data to analyze without having
   to prompt the user
   Remember to update varname_from_fname if you add a new file or if you change the name of the variable.

   fform i.e. the integer specification for data type is automatically initialized from varname.

SOURCE

672 subroutine fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg,ebands)
673 
674 !Arguments ------------------------------------
675 !scalars
676  integer,intent(in) :: iomode,cplex,nfft,nspden
677  character(len=*),intent(in) :: varname,path
678  type(hdr_type),intent(inout) :: hdr
679  type(crystal_t),intent(in) :: crystal
680  type(ebands_t),optional,intent(in) :: ebands
681  type(MPI_type),intent(in) :: mpi_enreg
682 !arrays
683  integer,intent(in) :: ngfft(18)
684  real(dp),intent(inout) :: datar(cplex*nfft,nspden)
685  !type(pawrhoij_type),optional,intent(inout) :: pawrhoij_all(hdr%usepaw*crystal%natom)
686 
687 !Local variables-------------------------------
688 !!scalars
689  integer,parameter :: master=0
690  integer :: n1,n2,n3,comm_fft,nproc_fft,me_fft,iarr,ierr,ispden,unt,mpierr,fform
691  integer :: i3_glob,my_iomode
692  integer(kind=XMPI_OFFSET_KIND) :: hdr_offset,my_offset,nfft_tot
693  integer :: ncid,ncerr
694  character(len=fnlen) :: file_etsf
695  real(dp) :: cputime,walltime,gflops
696  character(len=500) :: msg,errmsg
697  type(abifile_t) :: abifile
698 !arrays
699  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
700  integer(XMPI_OFFSET_KIND) :: bsize_frecord(nspden)
701 
702 ! *************************************************************************
703 
704  abifile = abifile_from_varname(varname)
705  if (abifile%fform == 0) then
706     ABI_ERROR(sjoin("Cannot find any abifile object associated to varname:", varname))
707  end if
708  ! Get fform from abifile. TODO: check file extension
709  fform = abifile%fform
710 
711  comm_fft = mpi_enreg%comm_fft; nproc_fft = xmpi_comm_size(comm_fft); me_fft = mpi_enreg%me_fft
712  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = n1*n2*n3
713 
714  ! Select iomode
715  ! Use Fortran IO if nproc_fft 1, in principle this is not needed because the
716  ! MPI-IO code should produce binary files that are readable with Fortran-IO
717  ! but it seems that NAG uses its own binary format
718  my_iomode = iomode
719  if (my_iomode /= IO_MODE_ETSF .and. nproc_fft == 1) my_iomode = IO_MODE_FORTRAN
720  if (nproc_fft > 1 .and. my_iomode == IO_MODE_FORTRAN) my_iomode = IO_MODE_MPI
721 
722  call wrtout(std_out, sjoin(ch10, "fftdatar_write: About to write data to:", path, "with iomode:",iomode2str(my_iomode)))
723  call cwtime(cputime, walltime, gflops, "start")
724 
725  ! Get MPI-FFT tables from input ngfft
726  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
727 
728  select case (my_iomode)
729  case (IO_MODE_FORTRAN)
730    ABI_CHECK(nproc_fft == 1, "MPI-IO must be enabled when FFT parallelism is used")
731    if (open_file(path, msg, newunit=unt, form='unformatted', status='unknown', action="write") /= 0) then
732      ABI_ERROR(msg)
733    end if
734    call hdr%fort_write(unt, fform, ierr)
735    ABI_CHECK(ierr==0,"ierr !=0")
736    do ispden=1,nspden
737      write(unt, err=10, iomsg=errmsg) (datar(iarr,ispden), iarr=1,cplex * nfft)
738    end do
739    close(unt, err=10, iomsg=errmsg)
740 
741    ! Write PAW rhoij
742    !call pawrhoij_io(hdr%pawrhoij,unit,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,headform,"Write")
743    !call pawrhoij_io(rhoij_ptr,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,&
744    !                 HDR_LATEST_HEADFORM,"Write",form="netcdf")
745 
746 #ifdef HAVE_MPI_IO
747  case (IO_MODE_MPI)
748    ! Find the first z-plane treated by this node.
749    ! WARNING: Here I assume that the z-planes in real space
750    ! are distributed in contiguous blocks (as usually done in MPI-FFT)
751    do i3_glob=1,n3
752      if (me_fft == fftn3_distrib(i3_glob)) exit
753    end do
754    ABI_CHECK(i3_glob /= n3 +1, "This processor does not have z-planes!")
755 
756    ! Master writes the header.
757    if (me_fft == master) call hdr%write_to_fname(path, fform)
758    call xmpi_barrier(comm_fft) ! TODO: Non-blocking barrier.
759 
760    call MPI_FILE_OPEN(comm_fft, path, MPI_MODE_RDWR, xmpio_info, unt, mpierr)
761    ABI_CHECK_MPI(mpierr,"MPI_FILE_OPEN")
762 
763    ! Skip the header and get the offset of the header
764    call hdr_mpio_skip(unt,fform,hdr_offset)
765    !write(std_out,*)"i3_glob, nfft, hdr_offset,",i3_glob,nfft,hdr_offset,fftn3_distrib == me_fft
766 
767    ! Each proc writes a contiguous slice of the nspden records.
768    ! my_offset is the position inside the Fortran record.
769    do ispden=1,nspden
770      my_offset = hdr_offset + xmpio_bsize_frm + ((ispden - 1) * 2 * xmpio_bsize_frm) + &
771      ((i3_glob-1) * cplex * n1 * n2 * xmpi_bsize_dp)  + ((ispden-1) * cplex * nfft_tot * xmpi_bsize_dp)
772      call MPI_FILE_WRITE_AT_ALL(unt,my_offset,datar(:,ispden),cplex*nfft,MPI_DOUBLE_PRECISION,MPI_STATUS_IGNORE,mpierr)
773      ABI_CHECK_MPI(mpierr,"MPI_FILE_WRITE_AT_ALL")
774    end do
775 
776    ! master writes the fortran record markers.
777    if (me_fft == master) then
778      bsize_frecord = cplex * nfft_tot * xmpi_bsize_dp
779 #if 1
780      my_offset = hdr_offset
781      do ispden=1,nspden
782        call xmpio_write_frm(unt,my_offset,xmpio_single,bsize_frecord(ispden),mpierr)
783        ABI_CHECK_MPI(mpierr,"xmpio_write_frm")
784      end do
785 #else
786      ! TODO: Understand why this code does not work!
787      call xmpio_write_frmarkers(unt,hdr_offset,xmpio_single,nspden,bsize_frecord,ierr)
788      ABI_CHECK(ierr==0, "xmpio_write_frmarkers")
789 #endif
790    end if
791 
792    call MPI_FILE_CLOSE(unt,mpierr)
793    ABI_CHECK_MPI(mpierr,"FILE_CLOSE!")
794 
795    ! Add full pawrhoij datastructure at the end of the file.
796    !if (present(pawrhoij_all) .and. me_fft == master .and. hdr%usepaw == 1) then
797    !  if (open_file(path, msg, newunit=unt, form='unformatted', status='old', action="write", access="append") /= 0) then
798    !    ABI_ERROR(msg)
799    !  end if
800    !  call pawrhoij_io(pawrhoij_all,un,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Write")
801    !  close(unt)
802    !end if
803 #endif
804 
805  case (IO_MODE_ETSF)
806    file_etsf = nctk_ncify(path)
807 
808    ! Write datar.
809    ncerr = nctk_write_datar(varname,file_etsf,ngfft,cplex,nfft,nspden, &
810    comm_fft,fftn3_distrib,ffti3_local,datar,action="create")
811    NCF_CHECK(ncerr)
812    call xmpi_barrier(comm_fft)
813 
814    ! Master writes the header.
815    if (xmpi_comm_rank(comm_fft) == master) then
816      NCF_CHECK(nctk_open_modify(ncid, file_etsf, xmpi_comm_self))
817      NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.))
818      ! Add information on the crystalline structure.
819      NCF_CHECK(crystal%ncwrite(ncid))
820      if (present(ebands)) then
821        NCF_CHECK(ebands_ncwrite(ebands, ncid))
822      end if
823 
824      ! Add full pawrhoij datastructure.
825      !if (present(pawrhoij_all) .and. me_fft == master .and. hdr%usepaw == 1) then
826      !  call pawrhoij_io(pawrhoij_all,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Write", form="netcdf")
827      !end if
828 
829      NCF_CHECK(nf90_close(ncid))
830    end if
831 
832  case default
833    ABI_ERROR(sjoin("Wrong iomode:",itoa(my_iomode)))
834  end select
835 
836  call cwtime_report(" IO operation", cputime, walltime, gflops)
837 
838  return
839 
840  ! Handle Fortran IO error
841 10 continue
842  ABI_ERROR(errmsg)
843 
844 end subroutine fftdatar_write

m_ioarr/fftdatar_write_from_hdr [ Functions ]

[ Top ] [ m_ioarr ] [ Functions ]

NAME

 fftdatar_write_from_hdr

FUNCTION

 Write an array in real space on the FFT box to file.
 crystal and ebands are constructed from the Abinit header.

TODO

 This routine will be removed when crystal_t and ebands_t will become standard objects
 available in the GS/DFPT part.

INPUTS

 [eigen](mband*hdr%nkpt*hdr%nsppol)=GS eigenvalues
 See fftdatar_write for the meaning of the other variables.

OUTPUT

SOURCE

869 subroutine fftdatar_write_from_hdr(varname,path,iomode,hdr,ngfft,cplex,nfft,nspden,datar,mpi_enreg,eigen)
870 
871 !Arguments ------------------------------------
872 !scalars
873  integer,intent(in) :: iomode,cplex,nfft,nspden
874  character(len=*),intent(in) :: varname,path
875  type(hdr_type),intent(inout) :: hdr
876  type(MPI_type),intent(in) :: mpi_enreg
877 !arrays
878  integer,intent(in) :: ngfft(18)
879  real(dp),intent(inout) :: datar(cplex*nfft,nspden)
880  real(dp),optional,intent(in) :: eigen(:)
881 
882 !Local variables-------------------------------
883 !!scalars
884  integer :: mband
885  type(crystal_t) :: crystal
886  type(ebands_t) :: ebands
887 !arrays
888  real(dp),allocatable :: ene3d(:,:,:)
889 
890 ! *************************************************************************
891 
892  crystal = hdr%get_crystal()
893 
894  if (present(eigen)) then
895      mband = maxval(hdr%nband)
896      ABI_CHECK(size(eigen) ==  mband * hdr%nkpt * hdr%nsppol, "Wrong size(eigen)")
897      ABI_MALLOC(ene3d, (mband, hdr%nkpt, hdr%nsppol))
898      call unpack_eneocc(hdr%nkpt, hdr%nsppol, mband, hdr%nband, eigen, ene3d)
899      ebands = ebands_from_hdr(hdr, mband, ene3d)
900      ABI_FREE(ene3d)
901 
902     call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg,ebands=ebands)
903     call ebands_free(ebands)
904  else
905     call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg)
906  end if
907 
908  call crystal%free()
909 
910 end subroutine fftdatar_write_from_hdr

m_ioarr/fort_denpot_skip [ Functions ]

[ Top ] [ m_ioarr ] [ Functions ]

NAME

  fort_denpot_skip

FUNCTION

  Skip the header and the DEN/POT records. Mainly used to append data to a pre-existent file.
  Return exit code.

INPUTS

  unit=Fortran unit number (already opened in the caller).
  msg=Error message if ierr /= 0

SOURCE

1251 integer function fort_denpot_skip(unit, msg) result(ierr)
1252 
1253 !Arguments ------------------------------------
1254  integer,intent(in) :: unit
1255  character(len=*),intent(out) :: msg
1256 
1257 !Local variables-------------------------------
1258  integer :: ii,fform,nspden
1259  type(hdr_type) :: hdr
1260 
1261 ! *********************************************************************
1262 
1263  ierr = 1
1264  call hdr_fort_read(hdr, unit, fform)
1265  if (fform == 0) then
1266     msg = "hdr_fort_read returned fform == 0"; return
1267  end if
1268 
1269  nspden = hdr%nspden
1270  call hdr%free()
1271 
1272  ! Skip the records with v1.
1273  do ii=1,nspden
1274    read(unit, iostat=ierr, iomsg=msg)
1275    if (ierr /= 0) return
1276  end do
1277 
1278  ierr = 0
1279 
1280 end function fort_denpot_skip

m_ioarr/ioarr [ Functions ]

[ Top ] [ m_ioarr ] [ Functions ]

NAME

 ioarr

FUNCTION

 Read or write rho(r) or v(r), either ground-state or response-functions.
 If ground-state, these arrays are real, if response-functions, these arrays are complex.
 (in general, an array stored in unformatted form on a real space fft grid).
 rdwr=1 to read, 2 to write

 This subroutine should be called only by one processor in the writing mode

INPUTS

 (some may be output)
 accessfil=
    0 for FORTRAN_IO
    3 for ETSF_IO
    4 for MPI_IO
 cplex=1 for real array, 2 for complex
 nfft=Number of FFT points treated by this node.
 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 dtset <type(dataset_type)>=all input variables for this dataset
 fform=integer specification for data type:
   2 for wf; 52 for density; 102 for potential
   old format (prior to ABINITv2.0): 1, 51 and 101.
 fildata=file name
 hdr <type(hdr_type)>=the header of wf, den and pot files
  if rdwr=1 , used to compare with the hdr of the read disk file
  if rdwr=2 , used as the header of the written disk file
 mpi_enreg=information about MPI parallelization
 rdwr=choice parameter, see above
 rdwrpaw=1 only if rhoij PAW quantities have to be read (if rdwr=1)
 [single_proc]=True if only ONE MPI process is calling this routine. This usually happens when
   master calls ioarr to read data that is then broadcasted in the caller. Default: False.
   Note that singleproc is not compatible with FFT parallelism because nfft is assumed to be
   the total number of points in the FFT mesh.

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output
 arr(cplex*nfft,nspden)=array on real space grid, returned for rdwr=1, input for rdwr=2
 etotal=total energy (Ha), returned for rdwr=1
 === if rdwrpaw/=0 ===
  pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data

SOURCE

128 subroutine ioarr(accessfil,arr,dtset,etotal,fform,fildata,hdr,mpi_enreg, &
129 &                ngfft,cplex,nfft,pawrhoij,rdwr,rdwrpaw,wvl_den,single_proc)
130 
131 !Arguments ------------------------------------
132 !scalars
133  integer,intent(in) :: accessfil,cplex,nfft,rdwr,rdwrpaw
134  integer,intent(inout) :: fform
135  real(dp),intent(inout) :: etotal
136  character(len=*),intent(in) :: fildata
137  logical,optional,intent(in) :: single_proc
138  type(MPI_type),intent(inout) :: mpi_enreg
139  type(dataset_type),intent(in) :: dtset
140  type(hdr_type),intent(inout) :: hdr
141  type(wvl_denspot_type),optional, intent(in) :: wvl_den
142 !arrays
143  integer,intent(in) :: ngfft(18)
144  real(dp),intent(inout),target :: arr(cplex*nfft,dtset%nspden)
145  type(pawrhoij_type),intent(inout) :: pawrhoij(:)
146 
147 !Local variables-------------------------------
148  integer :: ncid,ncerr
149  character(len=fnlen) :: file_etsf
150 #ifdef HAVE_BIGDFT
151  integer :: i,i1,i2,i3,ia,ind,n1,n2,n3
152  integer :: zindex,zstart,zstop
153 #endif
154 !scalars
155  integer,parameter :: master=0
156  logical,parameter :: ALLOW_FFTINTERP=.True.
157  logical :: need_fftinterp,icheck_fft,qeq0
158  integer :: in_unt,out_unt,nfftot_in,nfftot_out,nspden,ncplxfft
159  integer :: iomode,fform_dum,iarr,ierr,ispden,me,me_fft,comm_fft
160  integer :: comm_cell,usewvl,unt
161  integer :: restart,restartpaw,spaceComm,spaceComm_io
162  real(dp) :: cputime,walltime,gflops
163  character(len=500) :: msg,errmsg
164  character(len=fnlen) :: my_fildata
165  character(len=nctk_slen) :: varname
166  type(hdr_type) :: hdr0
167  type(wffile_type) :: wff
168  type(MPI_type) :: MPI_enreg_seq
169 !arrays
170  integer :: ngfft_in(18),ngfft_out(18)
171  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
172  real(dp), ABI_CONTIGUOUS pointer :: arr_file(:,:),my_density(:,:)
173  real(dp),allocatable :: rhor_file(:,:),rhog_in(:,:),rhor_out(:,:),rhog_out(:,:)
174 
175 ! *************************************************************************
176 
177  DBG_ENTER("COLL")
178 
179  ncplxfft = cplex*nfft
180 
181  restartpaw=0
182  my_fildata = fildata
183  nspden = dtset%nspden; usewvl = dtset%usewvl
184 
185  ! Check validity of arguments--only rho(r) (51,52) and V(r) (101,102) are presently supported
186  if ( (fform-1)/2 /=25 .and. (fform-1)/2 /=50 ) then
187    write(msg,'(a,i0,a)')' Input fform= ',fform,' not allowed.'
188    ABI_BUG(msg)
189  end if
190 
191  ! Print input fform
192  if ( (fform-1)/2==25 .and. rdwr==1) then
193    msg = ' ioarr: reading density data '
194  else if ( (fform-1)/2==25 .and. rdwr==2) then
195    msg = ' ioarr: writing density data'
196  else if ( (fform-1)/2==50 .and. rdwr==1) then
197    msg = ' ioarr: reading potential data'
198  else if ( (fform-1)/2==50 .and. rdwr==2) then
199    msg = ' ioarr: writing potential data'
200  end if
201  call wrtout(std_out,msg)
202 
203  call wrtout(std_out, 'ioarr: file name is: '//TRIM(fildata))
204 
205  if (accessfil == IO_MODE_ETSF) then ! Initialize filename in case of ETSF file.
206    file_etsf = nctk_ncify(fildata)
207    call wrtout(std_out,sjoin('file name for ETSF access: ', file_etsf))
208  end if
209 
210 !Some definitions for MPI-IO access
211  spaceComm = mpi_enreg%comm_cell
212  comm_cell = mpi_enreg%comm_cell
213  comm_fft = mpi_enreg%comm_fft
214 
215  if (accessfil == 4) then
216    iomode=IO_MODE_MPI
217    if (rdwr==1) then
218      spaceComm=mpi_enreg%comm_cell
219    else
220      spaceComm=mpi_enreg%comm_fft
221    end if
222    me=xmpi_comm_rank(spaceComm)
223    if (mpi_enreg%nproc_fft>1) then
224      me_fft=mpi_enreg%me_fft
225      spaceComm_io=mpi_enreg%comm_fft
226    else
227      me_fft=0
228      spaceComm_io=xmpi_comm_self
229    end if
230  end if
231  if (usewvl==1) then
232    spaceComm=mpi_enreg%comm_cell
233    me=xmpi_comm_rank(spaceComm)
234  end if
235 
236  ! Change communicators and ranks if we are calling ioarr with one single processor.
237  if (present(single_proc)) then
238    if (single_proc) then
239      spaceComm = xmpi_comm_self
240      spaceComm_io = xmpi_comm_self
241      ABI_CHECK(mpi_enreg%nproc_fft == 1, "single_proc cannot be used when nproc_fft > 1")
242      comm_cell = xmpi_comm_self
243      comm_fft = xmpi_comm_self
244      me = 0
245    end if
246  end if
247 
248 !=======================================
249 !Handle input from disk file
250 !=======================================
251 
252  call cwtime(cputime, walltime, gflops, "start")
253 
254  if (rdwr==1) then
255    if (accessfil == 0 .or. accessfil == 4) then
256 
257      ! Here master checks if the input rho(r) is given on a FFT mesh that quals
258      ! the one used in the run. If not, we perform a Fourier interpolation, we write the
259      ! interpolated rho(r) to a temporary file and we use this file to restart.
260      if (ALLOW_FFTINTERP .and. usewvl==0) then
261        need_fftinterp = .False.; icheck_fft = .True.
262        ! only master checks the FFT mesh if MPI-IO. All processors read ngfft if Fortran-IO
263        ! Note that, when Fortran-IO is used, we don't know if the routine is called
264        ! by a single processor or by all procs in comm_cell hence we cannot broadcast my_fildata
265        ! inside spaceComm as done if accessfil == 4
266        if (accessfil == 4) icheck_fft = (xmpi_comm_rank(spaceComm)==master)
267 
268        if (icheck_fft) then
269          if (open_file(fildata,msg,newunit=in_unt,form='unformatted',status='old') /= 0) then
270            ABI_ERROR(msg)
271          end if
272 
273          call hdr_io(fform_dum,hdr0,rdwr,in_unt)
274          need_fftinterp = (ANY(hdr%ngfft/=hdr0%ngfft) )
275          qeq0=(hdr%qptn(1)**2+hdr%qptn(2)**2+hdr%qptn(3)**2<1.d-14)
276          ! FIXME: SHould handle double-grid if PAW
277          nfftot_in = product(hdr0%ngfft(1:3))
278          nfftot_out = product(hdr%ngfft(1:3))
279 
280          if (need_fftinterp) then
281            write(msg, "(2a,2(a,3(i0,1x)))")&
282             "Will perform Fourier interpolation since in and out ngfft differ",ch10,&
283             "ngfft in file: ",hdr0%ngfft,", expected ngfft: ",hdr%ngfft
284            ABI_WARNING(msg)
285 
286            ! Read rho(r) from file, interpolate it, write data and change fildata
287            ABI_MALLOC(rhor_file, (cplex*nfftot_in, hdr0%nspden))
288            ABI_MALLOC(rhog_in, (2, nfftot_in))
289            ABI_MALLOC(rhor_out, (cplex*nfftot_out, hdr0%nspden))
290            ABI_MALLOC(rhog_out, (2, nfftot_out))
291 
292            do ispden=1,hdr0%nspden
293              read(in_unt, err=10, iomsg=errmsg) (rhor_file(iarr,ispden), iarr=1,cplex*nfftot_in)
294            end do
295 
296            ngfft_in = dtset%ngfft; ngfft_out = dtset%ngfft
297            ngfft_in(1:3) = hdr0%ngfft(1:3); ngfft_out(1:3) = hdr%ngfft(1:3)
298            ngfft_in(4:6) = hdr0%ngfft(1:3); ngfft_out(4:6) = hdr%ngfft(1:3)
299            ngfft_in(9:18) = 0; ngfft_out(9:18) = 0
300            ngfft_in(10) = 1; ngfft_out(10) = 1
301 
302            call initmpi_seq(MPI_enreg_seq)
303            ! Which one is coarse? Note that this part is not very robust and can fail!
304            if (ngfft_in(2) * ngfft_in(3) < ngfft_out(2) * ngfft_out(3)) then
305              call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft_in(2),ngfft_in(3),'all')
306              call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfft_out(2),ngfft_out(3),'all')
307            else
308              call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfft_in(2),ngfft_in(3),'all')
309              call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft_out(2),ngfft_out(3),'all')
310            end if
311 
312            call fourier_interpol(cplex,hdr0%nspden,0,0,nfftot_in,ngfft_in,nfftot_out,ngfft_out,&
313             MPI_enreg_seq,rhor_file,rhor_out,rhog_in,rhog_out)
314 
315            call destroy_mpi_enreg(MPI_enreg_seq)
316 
317            ! MG Hack: Change fildata so that we will use this file to read the correct rho(r)
318            ! FIXME: This should be done in a cleaner way!
319            my_fildata = trim(fildata)//"__fftinterp_rhor__"
320            if (my_fildata == fildata) my_fildata = "__fftinterp_rhor__"
321            if (open_file(my_fildata,msg,newunit=out_unt,form='unformatted',status='unknown') /= 0) then
322              ABI_ERROR(msg)
323            end if
324            call hdr_io(fform_dum,hdr,2,out_unt)
325            do ispden=1,hdr0%nspden
326              write(out_unt, err=10, iomsg=errmsg) (rhor_out(iarr,ispden),iarr=1,cplex*nfftot_out)
327            end do
328            close(out_unt)
329 
330            ABI_FREE(rhor_file)
331            ABI_FREE(rhog_in)
332            ABI_FREE(rhor_out)
333            ABI_FREE(rhog_out)
334          end if ! need_fftinterp
335 
336          call hdr0%free()
337          close(in_unt, err=10, iomsg=errmsg)
338        end if ! master
339        if (accessfil == 4) call xmpi_bcast(my_fildata,master,spaceComm,ierr)
340      end if
341 
342      if (accessfil == 4) then
343        unt = get_unit()
344        call WffOpen(iomode,spaceComm,my_fildata,ierr,wff,0,me,unt,spaceComm_io)
345        call hdr_io(fform_dum,hdr0,rdwr,wff)
346        ! Compare the internal header and the header from the file
347        call hdr_check(fform,fform_dum,hdr,hdr0,'COLL',restart,restartpaw)
348 
349      else
350        if (open_file(my_fildata, msg, newunit=unt, form="unformatted", status="old", action="read") /= 0) then
351          ABI_ERROR(msg)
352        end if
353        ! Initialize hdr0, thanks to reading of unwff1
354        call hdr_io(fform_dum,hdr0,rdwr,unt)
355        ! Compare the internal header and the header from the file
356        call hdr_check(fform,fform_dum,hdr,hdr0,'COLL',restart,restartpaw)
357      end if
358      etotal=hdr0%etot
359 
360      ! NOTE: should check that restart is possible !!
361      !call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
362 
363      ! If nspden[file] /= nspden, need a temporary array
364      if (hdr0%nspden/=nspden) then
365        ABI_MALLOC(arr_file,(cplex*nfft,hdr0%nspden))
366      else
367        arr_file => arr
368      end if
369 
370      ! Read data
371      do ispden=1,hdr0%nspden
372        if(accessfil == 4) then
373          call xderiveRRecInit(wff,ierr)
374          call xderiveRead(wff,arr(1:ncplxfft,ispden),ncplxfft,spaceComm_io,ierr)
375          call xderiveRRecEnd(wff,ierr)
376 
377          !call xmpio_read_dp(mpi_fh,offset,sc_mode,ncount,buf,fmarker,mpierr,advance)
378          !do idat=1,ndat
379          !  do i3=1,n3
380          !    if( fftn3_distrib(i3) == me_fft) then
381          !      i3_local = ffti3_local(i3)
382          !      i3_ldat = i3_local + (idat - 1) * nd3proc
383          !      do i2=1,n2
384          !        frbase=n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * nfft
385          !        do i1=1,n1
386          !          fofr(i1+frbase)=workr(1,i1,i2,i3_ldat)
387          !        end do
388          !      end do
389          !    end if
390          !  end do
391          !end do
392 
393        else
394          read(unt, err=10, iomsg=errmsg) (arr_file(iarr,ispden),iarr=1,ncplxfft)
395        end if
396      end do
397 
398      if (accessfil == 4) then
399        call wffclose(wff,ierr)
400      else
401        close (unit=unt, err=10, iomsg=errmsg)
402      end if
403 
404    else if (accessfil == 3) then
405 
406      ! Read the header and broadcast it in comm_cell
407      ! FIXME: Use xmpi_comm_self for the time-being because, in loper, ioarr
408      ! is called by me==0
409      call hdr_read_from_fname(hdr0, file_etsf, fform_dum, comm_cell)
410      !call hdr_read_from_fname(hdr0, file_etsf, fform_dum, xmpi_comm_self)
411      ABI_CHECK(fform_dum/=0, "hdr_read_from_fname returned fform 0")
412 
413      ! Compare the internal header and the header from the file
414      call hdr_check(fform, fform_dum, hdr, hdr0, 'COLL', restart, restartpaw)
415 
416      ! If nspden[file] /= nspden, need a temporary array
417      if (hdr0%nspden /= nspden) then
418        ABI_MALLOC(arr_file,(cplex*nfft,hdr0%nspden))
419      else
420        arr_file => arr
421      end if
422 
423      if (usewvl == 1) then
424        ! Read the array
425        if (fform==52) then ! density
426          varname = "density"
427        else if (fform==102) then ! all potential forms!!!!
428          varname = "exchange_correlation_potential"
429        end if
430 
431        ! Open the file
432        NCF_CHECK(nctk_open_read(ncid, file_etsf, xmpi_comm_self))
433        NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, varname), arr_file))
434        NCF_CHECK(nf90_close(ncid))
435      else
436        ! Get MPI-FFT tables from input ngfft
437        call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
438 
439        ! Get the name of the netcdf variable from the ABINIT extension and read data.
440        varname = varname_from_fname(file_etsf)
441        ncerr = nctk_read_datar(file_etsf,varname,ngfft,cplex,nfft,hdr0%nspden,comm_fft,fftn3_distrib,ffti3_local,arr)
442        NCF_CHECK(ncerr)
443      end if
444 
445    else
446      write(msg,'(a,i0,a)')'Bad value for accessfil', accessfil, ' on read '
447      ABI_BUG(msg)
448    end if
449 
450    call wrtout(std_out,sjoin("data read from disk file: ", fildata))
451 
452    etotal=hdr0%etot
453 
454    ! Possibly need to convert the potential/density spin components
455    if (hdr0%nspden/=nspden) then
456      call denpot_spin_convert(arr_file,hdr0%nspden,arr,nspden,fform)
457      ABI_FREE(arr_file)
458    end if
459 
460    ! Eventually copy (or distribute) PAW data
461    if (rdwrpaw==1.and.restartpaw/=0) then
462      if (size(hdr0%pawrhoij) /= size(pawrhoij)) then
463        call pawrhoij_copy(hdr0%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
464                           keep_nspden=.true.)
465      else
466        call pawrhoij_copy(hdr0%pawrhoij,pawrhoij,keep_nspden=.true.)
467      end if
468    end if
469 
470    if (accessfil == 0 .or. accessfil == 3 .or. accessfil == 4) call hdr0%free()
471 
472  ! =======================================
473  ! Set up for writing data
474  ! =======================================
475  else if (rdwr==2) then
476 
477 !  In the wavelet case (isolated boundary counditions), the
478 !  arr array has a buffer that we need to remove.
479    if (usewvl == 1) then
480 #ifdef HAVE_BIGDFT
481      zindex = wvl_den%denspot%dpbox%nscatterarr(me, 3)
482      if (wvl_den%denspot%rhod%geocode == 'F') then
483        n1 = (wvl_den%denspot%dpbox%ndims(1) - 31) / 2
484        n2 = (wvl_den%denspot%dpbox%ndims(2) - 31) / 2
485        n3 = (wvl_den%denspot%dpbox%ndims(3) - 31) / 2
486        zstart = max(15 - zindex, 0)
487        zstop  = wvl_den%denspot%dpbox%nscatterarr(me, 2) + &
488 &       wvl_den%denspot%dpbox%nscatterarr(me, 4) - &
489 &       max(zindex + wvl_den%denspot%dpbox%nscatterarr(me, 2) &
490 &       - 2 * n3 - 15, 0)
491      else
492        ABI_ERROR('ioarr: WVL not implemented yet.')
493      end if
494      if (zstop - zstart + 1 > 0) then
495 !      Our slab contains (zstop - zstart + 1) elements
496        ABI_MALLOC(my_density,((n1*2)*(n2*2)*(zstop-zstart),nspden))
497 !      We copy the data except the buffer to my_density
498        ind = 0
499 
500        do i3 = zstart, zstop - 1, 1
501          ia = (i3 - 1) * dtset%ngfft(1) * dtset%ngfft(2)
502          do i2 = 0, 2 * n2 - 1, 1
503            i = ia + (i2 + 14) * dtset%ngfft(1) + 14
504            do i1 = 0, 2 * n1 - 1, 1
505              i   = i + 1
506              ind = ind + 1
507              my_density(ind, :) = arr(i, :)
508            end do
509          end do
510        end do
511      else
512        nullify(my_density)
513      end if
514 #else
515      BIGDFT_NOTENABLED_ERROR()
516      if(.false. .and. present(wvl_den))then
517        write(std_out,*)' One should not be here'
518      endif
519 #endif
520    end if
521 
522    ! Make sure ngfft agrees with hdr%ngfft.
523    if (usewvl == 0) then
524      if (any(ngfft(:3) /= hdr%ngfft(:3))) then
525        write(msg,"(2(a,3(1x,i0)))")"input ngfft: ",ngfft(:3),"differs from  hdr%ngfft: ",hdr%ngfft(:3)
526        ABI_ERROR(msg)
527      end if
528    end if
529 
530    if (accessfil == 0 .or. accessfil == 4) then
531      if(accessfil == 4) then
532        unt = get_unit()
533        call WffOpen(iomode,spaceComm,fildata,ierr,wff,0,me,unt)
534        call hdr_io(fform,hdr,rdwr,wff)
535      else
536        if (open_file(fildata, msg, newunit=unt, form='unformatted', status='unknown', action="write") /= 0) then
537          ABI_ERROR(msg)
538        end if
539 
540        ! Write header
541        call hdr_io(fform,hdr,rdwr,unt)
542      end if
543 
544      ! Write actual data
545      do ispden=1,nspden
546        if(accessfil == 4) then
547          call xderiveWRecInit(wff,ierr,me_fft)
548          call xderiveWrite(wff,arr(1:ncplxfft,ispden),ncplxfft,spaceComm_io,ierr)
549          call xderiveWRecEnd(wff,ierr,me_fft)
550        else
551          if (usewvl == 0) then
552            write(unt, err=10, iomsg=errmsg) (arr(iarr,ispden),iarr=1,ncplxfft)
553          else
554            write(unt, err=10, iomsg=errmsg) (my_density(iarr,ispden),iarr=1,size(my_density, 1))
555          end if
556        end if
557      end do
558 
559      if(accessfil == 4) then
560        call WffClose(wff,ierr)
561      else
562        close(unt, err=10, iomsg=errmsg)
563      end if
564 
565    else if ( accessfil == 3 ) then
566 
567      ! Master in comm_fft creates the file and writes the header.
568      if (xmpi_comm_rank(comm_fft) == 0) then
569        call hdr%write_to_fname(file_etsf, fform)
570      end if
571      call xmpi_barrier(comm_fft)
572 
573      ! Write the array
574      if (usewvl == 0) then
575        ! Get MPI-FFT tables from input ngfft
576        call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
577 
578        varname = varname_from_fname(file_etsf)
579        ncerr = nctk_write_datar(varname,file_etsf,ngfft,cplex,nfft,nspden,comm_fft,fftn3_distrib,ffti3_local,arr)
580        NCF_CHECK(ncerr)
581      else
582        NCF_CHECK(nctk_open_modify(ncid, file_etsf, xmpi_comm_self))
583 
584        if (fform==52) then ! density
585          varname = "density"
586          if (usewvl == 0) then
587            NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), arr))
588          else
589            NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), my_density))
590          end if
591        else if (fform==102) then ! all potential forms!!!!
592          varname = "exchange_correlation_potential"
593          NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), arr))
594        end if
595 
596        NCF_CHECK(nf90_close(ncid))
597      end if
598 
599    else
600      write(msg,'(a,i0,a)')'Bad value for accessfil', accessfil, ' on write '
601      ABI_ERROR(msg)
602    end if
603 
604    if (usewvl == 1 .and. associated(my_density)) then
605      ABI_FREE(my_density)
606    end if
607 
608    call wrtout(std_out,sjoin(' Data written to disk file:', fildata))
609 
610  else
611    write(msg,'(a,i0,a)')'Called with rdwr = ',rdwr,' not allowed.'
612    ABI_BUG(msg)
613  end if
614 
615  call cwtime_report(" IO operation", cputime, walltime, gflops)
616 
617  DBG_EXIT("COLL")
618 
619  return
620 
621  ! Handle Fortran IO error
622 10 continue
623  ABI_ERROR(errmsg)
624 
625 end subroutine ioarr

m_ioarr/read_rhor [ Functions ]

[ Top ] [ m_ioarr ] [ Functions ]

NAME

 read_rhor

FUNCTION

  Read the DEN file with name fname reporting the density on the real FFT mesh
  specified through the input variable ngfft. If the FFT mesh asked in input and that found
  on file differ, the routine performs a FFT interpolation and renormalize the density so that it
  integrates to the correct number of electrons. The interpolation is done only for NC.
  For PAW, this is not possible because one should include the onsite contribution so this task
  is delegated to the caller.

INPUTS

 fname=Name of the file
 cplex=1 if array is real, 2 if complex e.g. DFPT density.
 nspden=Number of spin density components.
 nfft=Number of FFT points (treated by this processor)
 ngfft(18)=Info on the FFT mesh.
 pawread= 1 if pawrhoij should be read from file, 0 otherwise. Meaningful only if usepaw==1.
 mpi_enreg<MPI_type>=Information about MPI parallelization
 comm=MPI communicator. See notes
 [check_hdr] <type(hdr_type)>=Optional. Used to compare with the hdr read from disk file
   The routine will abort if restart cannot be performed.
 [allow_interp]=If True, the density read from file will be interpolated if the mesh differs from the one
   expected by the caller. This option is usually used in **self-consistent** calculations.
   If False (default), the code stops if the two meshes are different.

OUTPUT

 orhor(cplex*nfft,nspden)=The density on the real space mesh.
 ohdr=Abinit header read from file.
 pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data. only
   if pawread==1. The arrays is supposed to be already allocated in the caller and its
   size must be consistent with the MPI communicator comm.

NOTES

   if xmpi_comm_size(comm) == 1, nfft shall be equal to nfftot, and len(pawrhoij) == natom
   This means that one can call this routine with

     if (xmpi_comm_rank(comm) == 0) call read_rhor(...., comm=xmpi_comm_self)

   to get the full array and pawrhoij(natom) on the master node.

   if xmpi_comm_size(comm) > 1, nfft represents the number of FFT points treated by this processor,
   and pawrhoij is dimensioned with my_natom
   All the processors inside comm and comm_atom should call this routine.

SOURCE

 963 subroutine read_rhor(fname, cplex, nspden, nfft, ngfft, pawread, mpi_enreg, orhor, ohdr, pawrhoij, comm, &
 964                      check_hdr, allow_interp) ! Optional
 965 
 966 !Arguments ------------------------------------
 967 !scalars
 968  integer,intent(in) :: cplex,nfft,nspden,pawread,comm
 969  character(len=*),intent(in) :: fname
 970  type(MPI_type),intent(in) :: mpi_enreg
 971  type(hdr_type),intent(out) :: ohdr
 972  type(hdr_type),optional,intent(in) :: check_hdr
 973  logical,optional,intent(in) :: allow_interp
 974 !arrays
 975  integer,intent(in) :: ngfft(18)
 976  real(dp),intent(out) :: orhor(cplex*nfft,nspden)
 977  type(pawrhoij_type),intent(inout) :: pawrhoij(:)
 978 
 979 !Local variables-------------------------------
 980 !scalars
 981  integer,parameter :: master=0
 982  integer :: unt,fform,iomode,my_rank,mybase,globase,cplex_file
 983  integer :: ispden,ifft,nfftot_file,nprocs,ierr,i1,i2,i3,i3_local,n1,n2,n3
 984  integer,parameter :: fform_den=52
 985  integer :: restart, restartpaw
 986  integer :: ncerr
 987  real(dp) :: ratio,ucvol
 988  real(dp) :: cputime,walltime,gflops
 989  logical :: need_interp,have_mpifft,allow_interp__
 990  character(len=500) :: msg,errmsg
 991  character(len=fnlen) :: my_fname
 992  character(len=nctk_slen) :: varname
 993 !arrays
 994  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
 995  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2)
 996  real(dp),allocatable :: rhor_file(:,:),rhor_tmp(:,:)
 997  type(pawrhoij_type),allocatable :: pawrhoij_file(:)
 998 
 999 ! *************************************************************************
1000 
1001  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1002  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); have_mpifft = (nfft /= product(ngfft(1:3)))
1003  allow_interp__ = .False.; if (present(allow_interp)) allow_interp__ = allow_interp
1004 
1005  call timab(1280,1,tsec)
1006  call wrtout(std_out, sjoin(" About to read data(r) from:", fname), do_flush=.True.)
1007  call cwtime(cputime, walltime, gflops, "start")
1008 
1009  ! Master node opens the file, read the header and the FFT data
1010  ! This approach facilitates the interpolation of the density if in_ngfft(1:3) /= file_ngfft(1:3)
1011  if (my_rank == master) then
1012    my_fname = fname
1013    if (nctk_try_fort_or_ncfile(my_fname, msg) /= 0 ) then
1014      ABI_ERROR(msg)
1015    end if
1016 
1017    iomode = iomode_from_fname(my_fname)
1018    select case (iomode)
1019 
1020    case (IO_MODE_FORTRAN, IO_MODE_MPI)
1021      if (open_file(my_fname, msg, newunit=unt, form='unformatted', status='old', action="read") /= 0) then
1022        ABI_ERROR(msg)
1023      end if
1024 
1025      call hdr_fort_read(ohdr, unt, fform)
1026 
1027      ! Check important dimensions.
1028      ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname))
1029      !if (fform /= fform_den) then
1030      !  write(msg, "(3a, 2(a, i0))")' File: ',trim(my_fname),ch10,' is not a density file. fform: ',fform,", expecting: ", fform_den
1031      !  ABI_WARNING(msg)
1032      !end if
1033      cplex_file = 1
1034      if (ohdr%pertcase /= 0) then
1035        cplex_file = 2; if (ohdr%qptn(1)**2 + ohdr%qptn(2)**2 + ohdr%qptn(3)**2 <1.d-14) cplex_file= 1
1036      end if
1037      ABI_CHECK(cplex_file == cplex, "cplex_file != cplex")
1038 
1039      ! Read FFT array (full box)
1040      nfftot_file = product(ohdr%ngfft(:3))
1041      ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
1042      do ispden=1,ohdr%nspden
1043        read(unt, err=10, iomsg=errmsg) (rhor_file(ifft,ispden), ifft=1,cplex*nfftot_file)
1044      end do
1045      close(unt)
1046 
1047    case (IO_MODE_ETSF)
1048      NCF_CHECK(nctk_open_read(unt, my_fname, xmpi_comm_self))
1049      call hdr_ncread(ohdr, unt, fform)
1050 
1051      ! Check important dimensions.
1052      ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname))
1053      !if (fform /= fform_den) then
1054      !  write(msg, "(2a, 2(a, i0))")' File: ',trim(my_fname),' is not a density file: fform= ',fform,", expecting:", fform_den
1055      !  ABI_WARNING(msg)
1056      !end if
1057 
1058      cplex_file = 1
1059      if (ohdr%pertcase /= 0) then
1060        cplex_file = 2; if (ohdr%qptn(1)**2 + ohdr%qptn(2)**2 + ohdr%qptn(3)**2 <1.d-14) cplex_file= 1
1061      end if
1062      ABI_CHECK(cplex_file == cplex, "cplex_file != cplex")
1063 
1064      ! Read FFT array (full box)
1065      nfftot_file = product(ohdr%ngfft(:3))
1066      ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
1067 
1068      varname = varname_from_fname(my_fname)
1069      ncerr= nf90_get_var(unt, nctk_idname(unt, varname), rhor_file, &
1070                         count=[cplex, ohdr%ngfft(1), ohdr%ngfft(2), ohdr%ngfft(3), ohdr%nspden])
1071      NCF_CHECK(ncerr)
1072      NCF_CHECK(nf90_close(unt))
1073 
1074    case default
1075      ABI_ERROR(sjoin("Wrong iomode:", itoa(iomode)))
1076    end select
1077 
1078    need_interp = any(ohdr%ngfft(1:3) /= ngfft(1:3))
1079    if (need_interp) then
1080      msg = sjoin("Different FFT meshes. Caller expects:", ltoa(ngfft(1:3)), &
1081                  ". File: ", ltoa(ohdr%ngfft(1:3)), ". Need to perform interpolation.")
1082      ABI_COMMENT(msg)
1083      if (.not. allow_interp__) then
1084        write(msg, "(5a)") &
1085         " Cannot continue as allow_interp = .False. ", ch10, &
1086         " Please set ngfft to: ", trim(ltoa(ohdr%ngfft(1:3))), " in the input file"
1087        ABI_ERROR(msg)
1088      end if
1089 
1090      ABI_MALLOC(rhor_tmp, (cplex*product(ngfft(1:3)), ohdr%nspden))
1091      call timab(1281,1,tsec)
1092      call interpolate_denpot(cplex, ohdr%ngfft(1:3), ohdr%nspden, rhor_file, ngfft(1:3), rhor_tmp)
1093      call timab(1281,2,tsec)
1094 
1095      ohdr%ngfft(1:3) = ngfft(1:3)
1096      nfftot_file = product(ohdr%ngfft(:3))
1097      ABI_FREE(rhor_file)
1098      ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
1099      rhor_file = rhor_tmp
1100      ABI_FREE(rhor_tmp)
1101 
1102      ! Renormalize charge to avoid errors due to the interpolation.
1103      ! Do this only for NC since for PAW we should add the onsite contribution.
1104      ! This is left to the caller.
1105      !if (ohdr%usepaw == 0) then
1106      if (ohdr%usepaw == 0 .and. fform == fform_den) then
1107        call metric(gmet, gprimd, -1, rmet, ohdr%rprimd, ucvol)
1108        ratio = ohdr%nelect / (sum(rhor_file(:,1))*ucvol/ product(ngfft(1:3)))
1109        rhor_file = rhor_file * ratio
1110        write(msg,'(a,f8.2,a,f8.4)')' Expected nelect: ',ohdr%nelect,' renormalization ratio: ',ratio
1111        call wrtout(std_out,msg)
1112      end if
1113    end if ! need_interp
1114 
1115    ! Read PAW Rhoij
1116    if (ohdr%usepaw == 1) then
1117      ABI_MALLOC(pawrhoij_file, (ohdr%natom))
1118      call pawrhoij_nullify(pawrhoij_file)
1119      call pawrhoij_alloc(pawrhoij_file, ohdr%pawrhoij(1)%cplex_rhoij, ohdr%pawrhoij(1)%nspden, ohdr%pawrhoij(1)%nspinor, &
1120          ohdr%pawrhoij(1)%nsppol, ohdr%typat, lmnsize=ohdr%lmn_size, qphase=ohdr%pawrhoij(1)%qphase)
1121      call pawrhoij_copy(ohdr%pawrhoij, pawrhoij_file)
1122    end if
1123 
1124   ! Check that restart is possible !
1125   ! This check must be done here because we may have changed hdr% if need_interp
1126   if (present(check_hdr)) then
1127     ! FIXME: Temporary hack: fform_den to make hdr_check happy!
1128     call hdr_check(fform_den, fform_den, check_hdr, ohdr, "COLL", restart, restartpaw)
1129     !call hdr_check(fform_den, fform, check_hdr, ohdr, "COLL", restart, restartpaw)
1130   end if
1131 
1132  end if ! master
1133 
1134  if (nprocs == 1) then
1135    if (ohdr%nspden == nspden) then
1136      orhor = rhor_file
1137    else
1138      call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform)
1139    end if
1140    if (pawread == 1) call pawrhoij_copy(pawrhoij_file, pawrhoij, keep_nspden=.true.)
1141 
1142  else
1143    call ohdr%bcast(master, my_rank, comm)
1144    call xmpi_bcast(fform, master, comm, ierr)
1145 
1146    ! Eventually copy (or distribute) PAW data
1147    if (ohdr%usepaw == 1 .and. pawread == 1) then
1148      if (my_rank /= master) then
1149        ABI_MALLOC(pawrhoij_file, (ohdr%natom))
1150        call pawrhoij_nullify(pawrhoij_file)
1151        call pawrhoij_alloc(pawrhoij_file, ohdr%pawrhoij(1)%cplex_rhoij, ohdr%pawrhoij(1)%nspden, ohdr%pawrhoij(1)%nspinor, &
1152             ohdr%pawrhoij(1)%nsppol, ohdr%typat, lmnsize=ohdr%lmn_size, qphase=ohdr%pawrhoij(1)%qphase)
1153      end if
1154      if (size(ohdr%pawrhoij) /= size(pawrhoij)) then
1155        call pawrhoij_copy(ohdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
1156                           keep_nspden=.true.)
1157      else
1158        call pawrhoij_copy(ohdr%pawrhoij,pawrhoij, keep_nspden=.true.)
1159      end if
1160    end if
1161 
1162    if (my_rank /= master) then
1163      nfftot_file = product(ohdr%ngfft(1:3))
1164      ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
1165    end if
1166    call xmpi_bcast(rhor_file, master, comm,ierr)
1167 
1168    if (have_mpifft) then
1169      ! Extract slice treated by this MPI-FFT process.
1170      call ptabs_fourdp(mpi_enreg, ngfft(2), ngfft(3), fftn2_distrib, ffti2_local, fftn3_distrib, ffti3_local)
1171      if (ohdr%nspden==nspden) then
1172        do ispden=1,nspden
1173          do i3=1,n3
1174            if (fftn3_distrib(i3) /= mpi_enreg%me_fft) cycle
1175            i3_local = ffti3_local(i3)
1176            do i2=1,n2
1177              mybase = cplex * (n1 * (i2-1 + n2*(i3_local-1)))
1178              globase = cplex * (n1 * (i2-1 + n2*(i3-1)))
1179              do i1=1,n1*cplex
1180                orhor(i1+mybase,ispden) = rhor_file(i1+globase,ispden)
1181              end do
1182            end do
1183          end do
1184        end do
1185      else
1186        do i3=1,n3
1187          if (fftn3_distrib(i3) /= mpi_enreg%me_fft) cycle
1188          i3_local = ffti3_local(i3)
1189          do i2=1,n2
1190            mybase  = 1 + cplex * (n1 * (i2-1 + n2*(i3_local-1)))
1191            globase = 1 + cplex * (n1 * (i2-1 + n2*(i3-1)))
1192            call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform,&
1193                                     istart_in=globase,istart_out=mybase,nelem=n1*cplex)
1194          end do
1195        end do
1196      end if
1197    else
1198      if (ohdr%nspden==nspden) then
1199        orhor = rhor_file
1200      else
1201        call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform)
1202      end if
1203    end if
1204  end if ! nprocs > 1
1205 
1206  ABI_FREE(rhor_file)
1207 
1208  if (allocated(pawrhoij_file)) then
1209    call pawrhoij_free(pawrhoij_file)
1210    ABI_FREE(pawrhoij_file)
1211  end if
1212 
1213 ! Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities
1214 ! Add a small real to the magnetization
1215  if (nspden==4) orhor(:,4)=orhor(:,4)+tol14
1216  if (ohdr%usepaw==1.and.size(pawrhoij)>0) then
1217    if (pawrhoij(1)%nspden==4) then
1218      do i1=1,size(pawrhoij)
1219        pawrhoij(i1)%rhoijp(:,4)=pawrhoij(i1)%rhoijp(:,4)+tol10
1220      end do
1221    end if
1222  end if
1223 
1224  call timab(1280,2,tsec)
1225  call cwtime_report(" read_rhor", cputime, walltime, gflops)
1226  return
1227 
1228  ! Handle Fortran IO error
1229 10 continue
1230  ABI_ERROR(errmsg)
1231 
1232 end subroutine read_rhor