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-2018 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 .

PARENTS

SOURCE

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

PARENTS

      m_ioarr

CHILDREN

SOURCE

1428 subroutine denpot_spin_convert(denpot_in,nspden_in,denpot_out,nspden_out,fform,&
1429 &                              istart_in,istart_out,nelem) ! optional arguments
1430 
1431 
1432 !This section has been created automatically by the script Abilint (TD).
1433 !Do not modify the following lines by hand.
1434 #undef ABI_FUNC
1435 #define ABI_FUNC 'denpot_spin_convert'
1436 !End of the abilint section
1437 
1438  implicit none
1439 
1440 !Arguments ------------------------------------
1441 !scalars
1442  integer,intent(in) :: nspden_in,nspden_out,fform
1443  integer,intent(in),optional :: istart_in,istart_out,nelem
1444 !arrays
1445  real(dp),intent(in) :: denpot_in(:,:)
1446  real(dp),intent(out) :: denpot_out(:,:)
1447 
1448 !Local variables-------------------------------
1449  integer :: iend_in,iend_out,ispden,my_istart_in,my_istart_out,my_nelem
1450  character(len=500) :: msg
1451 
1452 ! *********************************************************************
1453 
1454 !Optional arguments
1455  my_istart_in=1;if (present(istart_in)) my_istart_in=istart_in
1456  my_istart_out=1;if (present(istart_out)) my_istart_out=istart_out
1457  iend_in=size(denpot_in,1) ; iend_out=size(denpot_out,1)
1458  my_nelem=min(iend_in-my_istart_in+1,iend_out-my_istart_out+1)
1459  if (present(nelem)) my_nelem=nelem
1460 
1461 !Checks
1462  if (size(denpot_in,2)/=nspden_in) then
1463    msg='size(denpot_in,2)/=nspden_in!'
1464    MSG_BUG(msg)
1465  end if
1466  if (size(denpot_out,2)/=nspden_out) then
1467    msg='size(denpot_out,2)/=nspden_out!'
1468    MSG_BUG(msg)
1469  end if
1470  if (my_istart_in+my_nelem-1>size(denpot_in,1)) then
1471    msg='istart_in+nelem>size(denpot_in,1)!'
1472    MSG_BUG(msg)
1473  end if
1474  if (my_istart_out+my_nelem-1>size(denpot_out,1)) then
1475    msg='istart_out+nelem>size(denpot_out,1)!'
1476    MSG_BUG(msg)
1477  end if
1478 
1479 !Simple copy if the number of spin-components is unchanged...
1480  if (nspden_in==nspden_out) then
1481    do ispden=1,nspden_in
1482      denpot_out(my_istart_out:my_istart_out+my_nelem-1,ispden)= &
1483 &      denpot_in(my_istart_in:my_istart_in+my_nelem-1,ispden)
1484    end do
1485    return
1486  end if
1487 
1488 !...otherwise, we need to convert.
1489  if ((fform-1)/2==25) then
1490 
1491 !  First case: DENSITY
1492 
1493    if      (nspden_in==1.and.nspden_out==2) then
1494      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1495      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
1496    else if (nspden_in==1.and.nspden_out==4) then
1497      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1498      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=zero
1499      denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1500      denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1501    else if (nspden_in==2.and.nspden_out==1) then
1502      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1503    else if (nspden_in==2.and.nspden_out==4) then
1504      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1505      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=zero
1506      denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1507      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 &
1508 &                                                        -denpot_in(my_istart_in:my_istart_in+my_nelem,1)
1509    else if (nspden_in==4.and.nspden_out==1) then
1510      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1511    else if (nspden_in==4.and.nspden_out==2) then
1512      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1513      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 &
1514 &                                                        +denpot_in(my_istart_in:my_istart_in+my_nelem-1,4)*half
1515    end if
1516 
1517  else
1518 
1519 !  Second case: POTENTIAL
1520 
1521    if      (nspden_in==1.and.nspden_out==2) then
1522      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1523      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1524    else if (nspden_in==1.and.nspden_out==4) then
1525      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1526      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1527      denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1528      denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1529    else if (nspden_in==2.and.nspden_out==1) then
1530      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 &
1531 &                                                        +denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*half
1532    else if (nspden_in==2.and.nspden_out==4) then
1533      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1534      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)
1535      denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1536      denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1537    else if (nspden_in==4.and.nspden_out==1) then
1538      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 &
1539 &                                                        +denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*half
1540    else if (nspden_in==4.and.nspden_out==2) then
1541      denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1542      denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)
1543    end if
1544 
1545  end if
1546 
1547 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.

PARENTS

      cut3d,m_ioarr,outscfcv,sigma

CHILDREN

SOURCE

708 subroutine fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg,ebands)
709 
710 
711 !This section has been created automatically by the script Abilint (TD).
712 !Do not modify the following lines by hand.
713 #undef ABI_FUNC
714 #define ABI_FUNC 'fftdatar_write'
715 !End of the abilint section
716 
717  implicit none
718 
719 !Arguments ------------------------------------
720 !scalars
721  integer,intent(in) :: iomode,cplex,nfft,nspden
722  character(len=*),intent(in) :: varname,path
723  type(hdr_type),intent(inout) :: hdr
724  type(crystal_t),intent(in) :: crystal
725  type(ebands_t),optional,intent(in) :: ebands
726  type(MPI_type),intent(in) :: mpi_enreg
727 !arrays
728  integer,intent(in) :: ngfft(18)
729  real(dp),intent(inout) :: datar(cplex*nfft,nspden)
730  !type(pawrhoij_type),optional,intent(inout) :: pawrhoij_all(hdr%usepaw*crystal%natom)
731 
732 !Local variables-------------------------------
733 !!scalars
734  integer,parameter :: master=0
735  integer :: n1,n2,n3,comm_fft,nproc_fft,me_fft,iarr,ierr,ispden,unt,mpierr,fform
736  integer :: i3_glob,my_iomode
737  integer(kind=XMPI_OFFSET_KIND) :: hdr_offset,my_offset,nfft_tot
738 #ifdef HAVE_NETCDF
739  integer :: ncid,ncerr
740  character(len=fnlen) :: file_etsf
741 #endif
742  real(dp) :: cputime,walltime,gflops
743  character(len=500) :: msg,errmsg
744  type(abifile_t) :: abifile
745 !arrays
746  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
747  integer(XMPI_OFFSET_KIND) :: bsize_frecord(nspden)
748 
749 ! *************************************************************************
750 
751  abifile = abifile_from_varname(varname)
752  if (abifile%fform == 0) then
753     MSG_ERROR(sjoin("Cannot find any abifile object associated to varname:", varname))
754  end if
755  ! Get fform from abifile. TODO: check file extension
756  fform = abifile%fform
757 
758  comm_fft = mpi_enreg%comm_fft; nproc_fft = xmpi_comm_size(comm_fft); me_fft = mpi_enreg%me_fft
759  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = n1*n2*n3
760 
761  ! Select iomode
762  ! Use Fortran IO if nproc_fft 1, in principle this is not needed because the
763  ! MPI-IO code should produce binary files that are readable with Fortran-IO
764  ! but it seems that NAG uses its own binary format
765  my_iomode = iomode
766  if (my_iomode /= IO_MODE_ETSF .and. nproc_fft == 1) my_iomode = IO_MODE_FORTRAN
767  if (nproc_fft > 1 .and. my_iomode == IO_MODE_FORTRAN) my_iomode = IO_MODE_MPI
768 
769  call wrtout(std_out, sjoin(" fftdatar_write: About to write data to:", path, "with iomode", iomode2str(my_iomode)))
770  call cwtime(cputime, walltime, gflops, "start")
771 
772  ! Get MPI-FFT tables from input ngfft
773  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
774 
775  select case (my_iomode)
776  case (IO_MODE_FORTRAN)
777    ABI_CHECK(nproc_fft == 1, "MPI-IO must be enabled when FFT parallelism is used")
778    if (open_file(path, msg, newunit=unt, form='unformatted', status='unknown', action="write") /= 0) then
779      MSG_ERROR(msg)
780    end if
781    call hdr_fort_write(hdr,unt,fform,ierr)
782    ABI_CHECK(ierr==0,"ierr !=0")
783    do ispden=1,nspden
784      write(unt, err=10, iomsg=errmsg) (datar(iarr,ispden), iarr=1,cplex * nfft)
785    end do
786    close(unt, err=10, iomsg=errmsg)
787 
788    ! Write PAW rhoij
789    !call pawrhoij_io(hdr%pawrhoij,unit,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,headform,"Write")
790    !call pawrhoij_io(rhoij_ptr,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,&
791    !                 HDR_LATEST_HEADFORM,"Write",form="netcdf")
792 
793 #ifdef HAVE_MPI_IO
794  case (IO_MODE_MPI)
795    ! Find the first z-plane treated by this node.
796    ! WARNING: Here I assume that the z-planes in real space
797    ! are distributed in contiguous blocks (as usually done in MPI-FFT)
798    do i3_glob=1,n3
799      if (me_fft == fftn3_distrib(i3_glob)) exit
800    end do
801    ABI_CHECK(i3_glob /= n3 +1, "This processor does not have z-planes!")
802 
803    ! Master writes the header.
804    if (me_fft == master) call hdr_write_to_fname(hdr,path,fform)
805    call xmpi_barrier(comm_fft) ! TODO: Non-blocking barrier.
806 
807    call MPI_FILE_OPEN(comm_fft, path, MPI_MODE_RDWR, xmpio_info, unt, mpierr)
808    ABI_CHECK_MPI(mpierr,"MPI_FILE_OPEN")
809 
810    ! Skip the header and get the offset of the header
811    call hdr_mpio_skip(unt,fform,hdr_offset)
812    !write(std_out,*)"i3_glob, nfft, hdr_offset,",i3_glob,nfft,hdr_offset,fftn3_distrib == me_fft
813 
814    ! Each proc writes a contiguous slice of the nspden records.
815    ! my_offset is the position inside the Fortran record.
816    do ispden=1,nspden
817      my_offset = hdr_offset + xmpio_bsize_frm + ((ispden - 1) * 2 * xmpio_bsize_frm) + &
818      ((i3_glob-1) * cplex * n1 * n2 * xmpi_bsize_dp)  + ((ispden-1) * cplex * nfft_tot * xmpi_bsize_dp)
819      call MPI_FILE_WRITE_AT_ALL(unt,my_offset,datar(:,ispden),cplex*nfft,MPI_DOUBLE_PRECISION,MPI_STATUS_IGNORE,mpierr)
820      ABI_CHECK_MPI(mpierr,"MPI_FILE_WRITE_AT_ALL")
821    end do
822 
823    ! master writes the fortran record markers.
824    if (me_fft == master) then
825      bsize_frecord = cplex * nfft_tot * xmpi_bsize_dp
826 #if 1
827      my_offset = hdr_offset
828      do ispden=1,nspden
829        call xmpio_write_frm(unt,my_offset,xmpio_single,bsize_frecord(ispden),mpierr)
830        ABI_CHECK_MPI(mpierr,"xmpio_write_frm")
831      end do
832 #else
833      ! TODO: Understand why this code does not work!
834      call xmpio_write_frmarkers(unt,hdr_offset,xmpio_single,nspden,bsize_frecord,ierr)
835      ABI_CHECK(ierr==0, "xmpio_write_frmarkers")
836 #endif
837    end if
838 
839    call MPI_FILE_CLOSE(unt,mpierr)
840    ABI_CHECK_MPI(mpierr,"FILE_CLOSE!")
841 
842    ! Add full pawrhoij datastructure at the end of the file.
843    !if (present(pawrhoij_all) .and. me_fft == master .and. hdr%usepaw == 1) then
844    !  if (open_file(path, msg, newunit=unt, form='unformatted', status='old', action="write", access="append") /= 0) then
845    !    MSG_ERROR(msg)
846    !  end if
847    !  call pawrhoij_io(pawrhoij_all,un,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Write")
848    !  close(unt)
849    !end if
850 #endif
851 
852 #ifdef HAVE_NETCDF
853  case (IO_MODE_ETSF)
854    file_etsf = nctk_ncify(path)
855 
856    ! Write datar.
857    ncerr = nctk_write_datar(varname,file_etsf,ngfft,cplex,nfft,nspden, &
858    comm_fft,fftn3_distrib,ffti3_local,datar,action="create")
859    NCF_CHECK(ncerr)
860    call xmpi_barrier(comm_fft)
861 
862    ! Master writes the header.
863    if (xmpi_comm_rank(comm_fft) == master) then
864      NCF_CHECK(nctk_open_modify(ncid, file_etsf, xmpi_comm_self))
865      NCF_CHECK(hdr_ncwrite(hdr, ncid, fform, nc_define=.True.))
866      ! Add information on the crystalline structure.
867      NCF_CHECK(crystal_ncwrite(crystal, ncid))
868      if (present(ebands)) then
869        NCF_CHECK(ebands_ncwrite(ebands, ncid))
870      end if
871 
872      ! Add full pawrhoij datastructure.
873      !if (present(pawrhoij_all) .and. me_fft == master .and. hdr%usepaw == 1) then
874      !  call pawrhoij_io(pawrhoij_all,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Write", form="netcdf")
875      !end if
876 
877      NCF_CHECK(nf90_close(ncid))
878    end if
879 #endif
880 
881  case default
882    MSG_ERROR(sjoin("Wrong iomode:",itoa(my_iomode)))
883  end select
884 
885  call cwtime(cputime, walltime, gflops, "stop")
886  write(msg,'(2(a,f9.1),a)')" IO operation completed. cpu_time: ",cputime," [s], walltime: ",walltime," [s]"
887  call wrtout(std_out, msg, "COLL", do_flush=.True.)
888 
889  return
890 
891  ! Handle Fortran IO error
892 10 continue
893  MSG_ERROR(errmsg)
894 
895 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

PARENTS

      dfpt_scfcv,scfcv

CHILDREN

SOURCE

925 subroutine fftdatar_write_from_hdr(varname,path,iomode,hdr,ngfft,cplex,nfft,nspden,datar,mpi_enreg,eigen)
926 
927 
928 !This section has been created automatically by the script Abilint (TD).
929 !Do not modify the following lines by hand.
930 #undef ABI_FUNC
931 #define ABI_FUNC 'fftdatar_write_from_hdr'
932 !End of the abilint section
933 
934  implicit none
935 
936 !Arguments ------------------------------------
937 !scalars
938  integer,intent(in) :: iomode,cplex,nfft,nspden
939  character(len=*),intent(in) :: varname,path
940  type(hdr_type),intent(inout) :: hdr
941  type(MPI_type),intent(in) :: mpi_enreg
942 !arrays
943  integer,intent(in) :: ngfft(18)
944  real(dp),intent(inout) :: datar(cplex*nfft,nspden)
945  real(dp),optional,intent(in) :: eigen(:)
946 
947 !Local variables-------------------------------
948 !!scalars
949  integer :: timrev,mband
950  type(crystal_t) :: crystal
951  type(ebands_t) :: ebands
952 !arrays
953  real(dp),allocatable :: ene3d(:,:,:)
954 
955 ! *************************************************************************
956 
957  timrev = 2; if (any(hdr%kptopt == [3, 4])) timrev = 1
958  call crystal_from_hdr(crystal, hdr, timrev)
959 
960  if (present(eigen)) then
961      mband = maxval(hdr%nband)
962      ABI_CHECK(size(eigen) ==  mband * hdr%nkpt * hdr%nsppol, "Wrong size(eigen)")
963      ABI_MALLOC(ene3d, (mband, hdr%nkpt, hdr%nsppol))
964      call unpack_eneocc(hdr%nkpt, hdr%nsppol, mband, hdr%nband, eigen, ene3d)
965      ebands = ebands_from_hdr(hdr, mband, ene3d)
966      ABI_FREE(ene3d)
967 
968     call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg,ebands=ebands)
969     call ebands_free(ebands)
970  else
971     call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg)
972  end if
973 
974  call crystal_free(crystal)
975 
976 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

PARENTS

CHILDREN

SOURCE

1354 integer function fort_denpot_skip(unit, msg) result(ierr)
1355 
1356 
1357 !This section has been created automatically by the script Abilint (TD).
1358 !Do not modify the following lines by hand.
1359 #undef ABI_FUNC
1360 #define ABI_FUNC 'fort_denpot_skip'
1361 !End of the abilint section
1362 
1363  implicit none
1364 
1365 !Arguments ------------------------------------
1366  integer,intent(in) :: unit
1367  character(len=*),intent(out) :: msg
1368 
1369 !Local variables-------------------------------
1370  integer :: ii,fform,nspden
1371  type(hdr_type) :: hdr
1372 
1373 ! *********************************************************************
1374 
1375  ierr = 1
1376  call hdr_fort_read(hdr, unit, fform)
1377  if (fform == 0) then
1378     msg = "hdr_fort_read returned fform == 0"; return
1379  end if
1380 
1381  nspden = hdr%nspden
1382  call hdr_free(hdr)
1383 
1384  ! Skip the records with v1.
1385  do ii=1,nspden
1386    read(unit, iostat=ierr, iomsg=msg)
1387    if (ierr /= 0) return
1388  end do
1389 
1390  ierr = 0
1391 
1392 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

PARENTS

      gstate,outscfcv

CHILDREN

SOURCE

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

INPUTS

 fname=Name of the density file
 cplex=1 if density 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.

PARENTS

      dfpt_looppert,dfptnl_loop,gstate,mrgscr,nonlinear,respfn,setup_positron
      sigma

CHILDREN

SOURCE

1034 subroutine read_rhor(fname, cplex, nspden, nfft, ngfft, pawread, mpi_enreg, orhor, ohdr, pawrhoij, comm, &
1035   check_hdr, allow_interp) ! Optional
1036 
1037 
1038 !This section has been created automatically by the script Abilint (TD).
1039 !Do not modify the following lines by hand.
1040 #undef ABI_FUNC
1041 #define ABI_FUNC 'read_rhor'
1042 !End of the abilint section
1043 
1044  implicit none
1045 
1046 !Arguments ------------------------------------
1047 !scalars
1048  integer,intent(in) :: cplex,nfft,nspden,pawread,comm
1049  character(len=*),intent(in) :: fname
1050  type(MPI_type),intent(in) :: mpi_enreg
1051  type(hdr_type),intent(out) :: ohdr
1052  type(hdr_type),optional,intent(in) :: check_hdr
1053  logical,optional,intent(in) :: allow_interp
1054 !arrays
1055  integer,intent(in) :: ngfft(18)
1056  real(dp),intent(out) :: orhor(cplex*nfft,nspden)
1057  type(pawrhoij_type),intent(inout) :: pawrhoij(:)
1058 
1059 !Local variables-------------------------------
1060 !scalars
1061  integer,parameter :: master=0,paral_kgb0=0
1062  integer :: unt,fform,iomode,my_rank,mybase,globase,cplex_file
1063 !integer :: optin,optout
1064  integer :: ispden,ifft,nfftot_file,nprocs,ierr,i1,i2,i3,i3_local,n1,n2,n3
1065  integer,parameter :: fform_den=52
1066  integer :: restart, restartpaw
1067 #ifdef HAVE_NETCDF
1068  integer :: ncerr
1069 #endif
1070  real(dp) :: ratio,ucvol
1071  real(dp) :: cputime,walltime,gflops
1072  logical :: need_interp,have_mpifft,allow_interp__
1073  character(len=500) :: msg,errmsg
1074  character(len=fnlen) :: my_fname
1075  character(len=nctk_slen) :: varname
1076  !type(mpi_type) :: mpi_enreg_seq
1077 !arrays
1078 !integer :: ngfft_file(18)
1079  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
1080  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1081 !real(dp) :: rhogdum(1)
1082  real(dp),allocatable :: rhor_file(:,:),rhor_tmp(:,:)
1083  type(pawrhoij_type),allocatable :: pawrhoij_file(:)
1084 
1085 ! *************************************************************************
1086 
1087  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1088  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); have_mpifft = (nfft /= product(ngfft(1:3)))
1089  allow_interp__ = .False.; if (present(allow_interp)) allow_interp__ = allow_interp
1090 
1091  call wrtout(std_out, sjoin(" About to read data(r) from:", fname), 'COLL')
1092  call cwtime(cputime, walltime, gflops, "start")
1093 
1094  ! Master node opens the file, read the header and the FFT data
1095  ! This approach facilitates the interpolation of the density if in_ngfft(1:3) /= file_ngfft(1:3)
1096  if (my_rank == master) then
1097    my_fname = fname
1098    if (nctk_try_fort_or_ncfile(my_fname, msg) /= 0 ) then
1099      MSG_ERROR(msg)
1100    end if
1101 
1102    iomode = iomode_from_fname(my_fname)
1103    select case (iomode)
1104 
1105    case (IO_MODE_FORTRAN, IO_MODE_MPI)
1106      if (open_file(my_fname, msg, newunit=unt, form='unformatted', status='old', action="read") /= 0) then
1107        MSG_ERROR(msg)
1108      end if
1109 
1110      call hdr_fort_read(ohdr, unt, fform)
1111 
1112      ! Check important dimensions.
1113      ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname))
1114      if (fform /= fform_den) then
1115        write(msg, "(2a, 2(a, i0))")' File: ',trim(my_fname),' is not a density file: fform= ',fform,", expecting:", fform_den
1116        MSG_WARNING(msg)
1117      end if
1118      cplex_file = 1
1119      if (ohdr%pertcase /= 0) then
1120        cplex_file = 2; if (ohdr%qptn(1)**2 + ohdr%qptn(2)**2 + ohdr%qptn(3)**2 <1.d-14) cplex_file= 1
1121      end if
1122      ABI_CHECK(cplex_file == cplex, "cplex_file != cplex")
1123 
1124      ! Read FFT array (full box)
1125      nfftot_file = product(ohdr%ngfft(:3))
1126      ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
1127      do ispden=1,ohdr%nspden
1128        read(unt, err=10, iomsg=errmsg) (rhor_file(ifft,ispden), ifft=1,cplex*nfftot_file)
1129      end do
1130      close(unt)
1131 
1132 #ifdef HAVE_NETCDF
1133    case (IO_MODE_ETSF)
1134      NCF_CHECK(nctk_open_read(unt, my_fname, xmpi_comm_self))
1135      call hdr_ncread(ohdr, unt, fform)
1136 
1137      ! Check important dimensions.
1138      ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname))
1139      if (fform /= fform_den) then
1140        write(msg, "(2a, 2(a, i0))")' File: ',trim(my_fname),' is not a density file: fform= ',fform,", expecting:", fform_den
1141        MSG_WARNING(msg)
1142      end if
1143 
1144      cplex_file = 1
1145      if (ohdr%pertcase /= 0) then
1146        cplex_file = 2; if (ohdr%qptn(1)**2 + ohdr%qptn(2)**2 + ohdr%qptn(3)**2 <1.d-14) cplex_file= 1
1147      end if
1148      ABI_CHECK(cplex_file == cplex, "cplex_file != cplex")
1149 
1150      ! Read FFT array (full box)
1151      nfftot_file = product(ohdr%ngfft(:3))
1152      ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
1153 
1154      varname = varname_from_fname(my_fname)
1155      ncerr= nf90_get_var(unt, nctk_idname(unt, varname), rhor_file, &
1156                         count=[cplex, ohdr%ngfft(1), ohdr%ngfft(2), ohdr%ngfft(3), ohdr%nspden])
1157      NCF_CHECK(ncerr)
1158      NCF_CHECK(nf90_close(unt))
1159 #endif
1160    case default
1161      MSG_ERROR(sjoin("Wrong iomode:", itoa(iomode)))
1162    end select
1163 
1164    need_interp = any(ohdr%ngfft(1:3) /= ngfft(1:3))
1165    if (need_interp .and. allow_interp__) then
1166      MSG_COMMENT("Real space meshes (DEN file, input rhor) are different. Will interpolate rhor(r).")
1167 
1168 #if 0
1169      ABI_CHECK(cplex == 1, "cplex != 1 not coded!")
1170      ngfft_file(1:3) = ohdr%ngfft(1:3)
1171      ngfft_file(4) = 2*(ngfft_file(1)/2)+1 ! 4:18 are used in fourdp
1172      ngfft_file(5) = 2*(ngfft_file(2)/2)+1
1173      ngfft_file(6) = ngfft_file(3)
1174      ngfft_file(7:18) = ngfft(7:18)
1175      optin  = 0 ! Input is taken from rhor
1176      optout = 0 ! Output is only in real space
1177 
1178      ! Fake MPI_type for the sequential part.
1179      !call initmpi_seq(MPI_enreg_seq)
1180      !call init_distribfft_seq(MPI_enreg_seq%distribfft, 'c', ngfftc(2), ngfftc(3), 'all')
1181      !call init_distribfft_seq(MPI_enreg_seq%distribfft, 'f', ngfftf(2), ngfftf(3), 'all')
1182 
1183      call fourier_interpol(cplex,ohdr%nspden,optin,optout,nfftot_file,ngfft_file,nfft,ngfft,&
1184        paral_kgb0,mpi_enreg,rhor_file,orhor,rhogdum,rhogdum)
1185 
1186      !call destroy_mpi_enreg(MPI_enreg_seq)
1187 #endif
1188 
1189      ABI_MALLOC(rhor_tmp, (cplex*product(ngfft(1:3)), ohdr%nspden))
1190      call interpolate_denpot(cplex, ohdr%ngfft(1:3), ohdr%nspden, rhor_file, ngfft(1:3), rhor_tmp)
1191 
1192      ohdr%ngfft(1:3) = ngfft(1:3)
1193      nfftot_file = product(ohdr%ngfft(:3))
1194      ABI_FREE(rhor_file)
1195      ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
1196      rhor_file = rhor_tmp
1197      ABI_FREE(rhor_tmp)
1198 
1199      ! Renormalize charge to avoid errors due to the interpolation.
1200      ! Do this only for NC since for PAW we should add the onsite contribution.
1201      ! This is left to the caller.
1202      if (ohdr%usepaw == 0) then
1203        call metric(gmet, gprimd, -1, rmet, ohdr%rprimd, ucvol)
1204        ratio = ohdr%nelect / (sum(rhor_file(:,1))*ucvol/ product(ngfft(1:3)))
1205        rhor_file = rhor_file * ratio
1206        write(msg,'(a,f8.2,a,f8.4)')' Expected nelect: ',ohdr%nelect,' renormalization ratio: ',ratio
1207        call wrtout(std_out,msg,'COLL')
1208      end if
1209    end if ! need_interp
1210 
1211    ! Read PAW Rhoij
1212    if (ohdr%usepaw == 1) then
1213      ABI_DT_MALLOC(pawrhoij_file, (ohdr%natom))
1214      call pawrhoij_nullify(pawrhoij_file)
1215      call pawrhoij_alloc(pawrhoij_file, ohdr%pawrhoij(1)%cplex, ohdr%pawrhoij(1)%nspden, ohdr%pawrhoij(1)%nspinor, &
1216          ohdr%pawrhoij(1)%nsppol, ohdr%typat, lmnsize=ohdr%lmn_size)
1217      call pawrhoij_copy(ohdr%pawrhoij, pawrhoij_file)
1218    end if
1219 
1220   ! Check that restart is possible !
1221   ! This check must be done here because we may have changed hdr% if need_interp
1222   if (present(check_hdr)) then
1223     ! FIXME: Temporary hack: fform_den to make hdr_check happy!
1224     call hdr_check(fform_den, fform_den, check_hdr, ohdr, "COLL", restart, restartpaw)
1225     !call hdr_check(fform_den, fform, check_hdr, ohdr, "COLL", restart, restartpaw)
1226   end if
1227 
1228 
1229  end if ! master
1230 
1231  if (nprocs == 1) then
1232    if (ohdr%nspden==nspden) then
1233      orhor = rhor_file
1234    else
1235      call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform)
1236    end if
1237    if (pawread == 1) call pawrhoij_copy(pawrhoij_file, pawrhoij, keep_nspden=.true.)
1238 
1239  else
1240    call hdr_bcast(ohdr, master, my_rank, comm)
1241    call xmpi_bcast(fform, master, comm, ierr)
1242 
1243    ! Eventually copy (or distribute) PAW data
1244    if (ohdr%usepaw == 1 .and. pawread == 1) then
1245      if (my_rank /= master) then
1246        ABI_DT_MALLOC(pawrhoij_file, (ohdr%natom))
1247        call pawrhoij_nullify(pawrhoij_file)
1248        call pawrhoij_alloc(pawrhoij_file, ohdr%pawrhoij(1)%cplex, ohdr%pawrhoij(1)%nspden, ohdr%pawrhoij(1)%nspinor, &
1249 &           ohdr%pawrhoij(1)%nsppol, ohdr%typat, lmnsize=ohdr%lmn_size)
1250      end if
1251      if (size(ohdr%pawrhoij) /= size(pawrhoij)) then
1252        call pawrhoij_copy(ohdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
1253 &                         keep_nspden=.true.)
1254      else
1255        call pawrhoij_copy(ohdr%pawrhoij,pawrhoij, keep_nspden=.true.)
1256      end if
1257    end if
1258 
1259    if (my_rank /= master) then
1260      nfftot_file = product(ohdr%ngfft(1:3))
1261      ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
1262    end if
1263    call xmpi_bcast(rhor_file, master, comm,ierr)
1264 
1265    if (have_mpifft) then
1266      ! Extract slice treated by this MPI-FFT process.
1267      call ptabs_fourdp(mpi_enreg, ngfft(2), ngfft(3), fftn2_distrib, ffti2_local, fftn3_distrib, ffti3_local)
1268      if (ohdr%nspden==nspden) then
1269        do ispden=1,nspden
1270          do i3=1,n3
1271            if (fftn3_distrib(i3) /= mpi_enreg%me_fft) cycle
1272            i3_local = ffti3_local(i3)
1273            do i2=1,n2
1274              mybase = cplex * (n1 * (i2-1 + n2*(i3_local-1)))
1275              globase = cplex * (n1 * (i2-1 + n2*(i3-1)))
1276              do i1=1,n1*cplex
1277                orhor(i1+mybase,ispden) = rhor_file(i1+globase,ispden)
1278              end do
1279            end do
1280          end do
1281        end do
1282      else
1283        do i3=1,n3
1284          if (fftn3_distrib(i3) /= mpi_enreg%me_fft) cycle
1285          i3_local = ffti3_local(i3)
1286          do i2=1,n2
1287            mybase  = 1 + cplex * (n1 * (i2-1 + n2*(i3_local-1)))
1288            globase = 1 + cplex * (n1 * (i2-1 + n2*(i3-1)))
1289            call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform,&
1290 &                  istart_in=globase,istart_out=mybase,nelem=n1*cplex)
1291          end do
1292        end do
1293      end if
1294    else
1295      if (ohdr%nspden==nspden) then
1296        orhor = rhor_file
1297      else
1298        call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform)
1299      end if
1300    end if
1301  end if ! nprocs > 1
1302 
1303  ABI_FREE(rhor_file)
1304 
1305  if (allocated(pawrhoij_file)) then
1306    call pawrhoij_free(pawrhoij_file)
1307    ABI_DT_FREE(pawrhoij_file)
1308  end if
1309 
1310 !Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities
1311 !  Add a small real to the magnetization
1312  if (nspden==4) orhor(:,4)=orhor(:,4)+tol14
1313  if (ohdr%usepaw==1.and.size(pawrhoij)>0) then
1314    if (pawrhoij(1)%nspden==4) then
1315      do i1=1,size(pawrhoij)
1316        pawrhoij(i1)%rhoijp(:,4)=pawrhoij(i1)%rhoijp(:,4)+tol10
1317      end do
1318    end if
1319  end if
1320 
1321  call cwtime(cputime, walltime, gflops, "stop")
1322  write(msg,'(2(a,f9.1),a)')" IO operation completed. cpu_time: ",cputime," [s], walltime: ",walltime," [s]"
1323  call wrtout(std_out, msg, "COLL", do_flush=.True.)
1324 
1325  return
1326 
1327  ! Handle Fortran IO error
1328 10 continue
1329  MSG_ERROR(errmsg)
1330 
1331 end subroutine read_rhor