TABLE OF CONTENTS
- ABINIT/m_ioarr
- m_ioarr/denpot_spin_convert
- m_ioarr/fftdatar_write
- m_ioarr/fftdatar_write_from_hdr
- m_ioarr/fort_denpot_skip
- m_ioarr/ioarr
- m_ioarr/read_rhor
ABINIT/m_ioarr [ 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-2024 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
1313 subroutine denpot_spin_convert(denpot_in,nspden_in,denpot_out,nspden_out,fform,& 1314 & istart_in,istart_out,nelem) ! optional arguments 1315 1316 !Arguments ------------------------------------ 1317 !scalars 1318 integer,intent(in) :: nspden_in,nspden_out,fform 1319 integer,intent(in),optional :: istart_in,istart_out,nelem 1320 !arrays 1321 real(dp),intent(in) :: denpot_in(:,:) 1322 real(dp),intent(out) :: denpot_out(:,:) 1323 1324 !Local variables------------------------------- 1325 integer :: iend_in,iend_out,ispden,my_istart_in,my_istart_out,my_nelem 1326 character(len=500) :: msg 1327 1328 ! ********************************************************************* 1329 1330 !Optional arguments 1331 my_istart_in=1;if (present(istart_in)) my_istart_in=istart_in 1332 my_istart_out=1;if (present(istart_out)) my_istart_out=istart_out 1333 iend_in=size(denpot_in,1) ; iend_out=size(denpot_out,1) 1334 my_nelem=min(iend_in-my_istart_in+1,iend_out-my_istart_out+1) 1335 if (present(nelem)) my_nelem=nelem 1336 1337 !Checks 1338 if (size(denpot_in,2)/=nspden_in) then 1339 msg='size(denpot_in,2)/=nspden_in!' 1340 ABI_BUG(msg) 1341 end if 1342 if (size(denpot_out,2)/=nspden_out) then 1343 msg='size(denpot_out,2)/=nspden_out!' 1344 ABI_BUG(msg) 1345 end if 1346 if (my_istart_in+my_nelem-1>size(denpot_in,1)) then 1347 msg='istart_in+nelem>size(denpot_in,1)!' 1348 ABI_BUG(msg) 1349 end if 1350 if (my_istart_out+my_nelem-1>size(denpot_out,1)) then 1351 msg='istart_out+nelem>size(denpot_out,1)!' 1352 ABI_BUG(msg) 1353 end if 1354 1355 !Simple copy if the number of spin-components is unchanged... 1356 if (nspden_in==nspden_out) then 1357 do ispden=1,nspden_in 1358 denpot_out(my_istart_out:my_istart_out+my_nelem-1,ispden)= & 1359 & denpot_in(my_istart_in:my_istart_in+my_nelem-1,ispden) 1360 end do 1361 return 1362 end if 1363 1364 !...otherwise, we need to convert. 1365 if ((fform-1)/2==25) then 1366 1367 ! First case: DENSITY 1368 1369 if (nspden_in==1.and.nspden_out==2) then 1370 denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1371 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 1372 else if (nspden_in==1.and.nspden_out==4) then 1373 denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1374 denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=zero 1375 denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero 1376 denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero 1377 else if (nspden_in==2.and.nspden_out==1) 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 else if (nspden_in==2.and.nspden_out==4) then 1380 denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1381 denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=zero 1382 denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero 1383 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 & 1384 & -denpot_in(my_istart_in:my_istart_in+my_nelem,1) 1385 else if (nspden_in==4.and.nspden_out==1) 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 else if (nspden_in==4.and.nspden_out==2) then 1388 denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1389 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 & 1390 & +denpot_in(my_istart_in:my_istart_in+my_nelem-1,4)*half 1391 end if 1392 1393 else 1394 1395 ! Second case: POTENTIAL 1396 1397 if (nspden_in==1.and.nspden_out==2) then 1398 denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1399 denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1400 else if (nspden_in==1.and.nspden_out==4) then 1401 denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1402 denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1403 denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero 1404 denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero 1405 else if (nspden_in==2.and.nspden_out==1) then 1406 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 & 1407 & +denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*half 1408 else if (nspden_in==2.and.nspden_out==4) then 1409 denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1410 denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2) 1411 denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero 1412 denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero 1413 else if (nspden_in==4.and.nspden_out==1) then 1414 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 & 1415 & +denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*half 1416 else if (nspden_in==4.and.nspden_out==2) then 1417 denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1) 1418 denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2) 1419 end if 1420 1421 end if 1422 1423 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
674 subroutine fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg,ebands) 675 676 !Arguments ------------------------------------ 677 !scalars 678 integer,intent(in) :: iomode,cplex,nfft,nspden 679 character(len=*),intent(in) :: varname,path 680 type(hdr_type),intent(inout) :: hdr 681 type(crystal_t),intent(in) :: crystal 682 type(ebands_t),optional,intent(in) :: ebands 683 type(MPI_type),intent(in) :: mpi_enreg 684 !arrays 685 integer,intent(in) :: ngfft(18) 686 real(dp),intent(inout) :: datar(cplex*nfft,nspden) 687 !type(pawrhoij_type),optional,intent(inout) :: pawrhoij_all(hdr%usepaw*crystal%natom) 688 689 !Local variables------------------------------- 690 !!scalars 691 integer,parameter :: master=0 692 integer :: n1,n2,n3,comm_fft,nproc_fft,me_fft,iarr,ierr,ii,ispden,unt,mpierr,fform 693 integer :: i3_glob,my_iomode 694 integer(kind=XMPI_OFFSET_KIND) :: hdr_offset,my_offset,nfft_tot 695 integer :: ncid,ncerr 696 character(len=fnlen) :: file_etsf 697 real(dp) :: cputime,walltime,gflops 698 character(len=500) :: msg,errmsg 699 type(abifile_t) :: abifile 700 !arrays 701 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:) 702 integer(XMPI_OFFSET_KIND) :: bsize_frecord(nspden) 703 704 ! ************************************************************************* 705 706 abifile = abifile_from_varname(varname) 707 if (abifile%fform == 0) then 708 ABI_ERROR(sjoin("Cannot find any abifile object associated to varname:", varname)) 709 end if 710 ! Get fform from abifile. TODO: check file extension 711 fform = abifile%fform 712 713 comm_fft = mpi_enreg%comm_fft; nproc_fft = xmpi_comm_size(comm_fft); me_fft = mpi_enreg%me_fft 714 n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = n1*n2*n3 715 716 ! Select iomode 717 ! Use Fortran IO if nproc_fft 1, in principle this is not needed because the 718 ! MPI-IO code should produce binary files that are readable with Fortran-IO 719 ! but it seems that NAG uses its own binary format 720 my_iomode = iomode 721 if (my_iomode /= IO_MODE_ETSF .and. nproc_fft == 1) my_iomode = IO_MODE_FORTRAN 722 if (nproc_fft > 1 .and. my_iomode == IO_MODE_FORTRAN) my_iomode = IO_MODE_MPI 723 724 call wrtout(std_out, sjoin(ch10, "fftdatar_write: About to write data to:", path, "with iomode:",iomode2str(my_iomode))) 725 call cwtime(cputime, walltime, gflops, "start") 726 727 ! Get MPI-FFT tables from input ngfft 728 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 729 730 select case (my_iomode) 731 case (IO_MODE_FORTRAN) 732 ABI_CHECK(nproc_fft == 1, "MPI-IO must be enabled when FFT parallelism is used") 733 if (open_file(path, msg, newunit=unt, form='unformatted', status='unknown', action="write") /= 0) then 734 ABI_ERROR(msg) 735 end if 736 call hdr%fort_write(unt, fform, ierr) 737 ABI_CHECK(ierr==0,"ierr !=0") 738 do ii=1,nspden 739 write(unt, err=10, iomsg=errmsg) (datar(iarr,ii), iarr=1,cplex * nfft) 740 end do 741 close(unt, err=10, iomsg=errmsg) 742 743 ! Write PAW rhoij 744 !call pawrhoij_io(hdr%pawrhoij,unit,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,headform,"Write") 745 !call pawrhoij_io(rhoij_ptr,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,& 746 ! HDR_LATEST_HEADFORM,"Write",form="netcdf") 747 748 #ifdef HAVE_MPI_IO 749 case (IO_MODE_MPI) 750 ! Find the first z-plane treated by this node. 751 ! WARNING: Here I assume that the z-planes in real space 752 ! are distributed in contiguous blocks (as usually done in MPI-FFT) 753 do i3_glob=1,n3 754 if (me_fft == fftn3_distrib(i3_glob)) exit 755 end do 756 ABI_CHECK(i3_glob /= n3 +1, "This processor does not have z-planes!") 757 758 ! Master writes the header. 759 if (me_fft == master) call hdr%write_to_fname(path, fform) 760 call xmpi_barrier(comm_fft) ! TODO: Non-blocking barrier. 761 762 call MPI_FILE_OPEN(comm_fft, path, MPI_MODE_RDWR, xmpio_info, unt, mpierr) 763 ABI_CHECK_MPI(mpierr,"MPI_FILE_OPEN") 764 765 ! Skip the header and get the offset of the header 766 call hdr_mpio_skip(unt,fform,hdr_offset) 767 !write(std_out,*)"i3_glob, nfft, hdr_offset,",i3_glob,nfft,hdr_offset,fftn3_distrib == me_fft 768 769 ! Each proc writes a contiguous slice of the nspden records. 770 ! my_offset is the position inside the Fortran record. 771 do ispden=1,nspden 772 my_offset = hdr_offset + xmpio_bsize_frm + ((ispden - 1) * 2 * xmpio_bsize_frm) + & 773 ((i3_glob-1) * cplex * n1 * n2 * xmpi_bsize_dp) + ((ispden-1) * cplex * nfft_tot * xmpi_bsize_dp) 774 call MPI_FILE_WRITE_AT_ALL(unt,my_offset,datar(:,ispden),cplex*nfft,MPI_DOUBLE_PRECISION,MPI_STATUS_IGNORE,mpierr) 775 ABI_CHECK_MPI(mpierr,"MPI_FILE_WRITE_AT_ALL") 776 end do 777 778 ! master writes the fortran record markers. 779 if (me_fft == master) then 780 bsize_frecord = cplex * nfft_tot * xmpi_bsize_dp 781 #if 1 782 my_offset = hdr_offset 783 do ispden=1,nspden 784 call xmpio_write_frm(unt,my_offset,xmpio_single,bsize_frecord(ispden),mpierr) 785 ABI_CHECK_MPI(mpierr,"xmpio_write_frm") 786 end do 787 #else 788 ! TODO: Understand why this code does not work! 789 call xmpio_write_frmarkers(unt,hdr_offset,xmpio_single,nspden,bsize_frecord,ierr) 790 ABI_CHECK(ierr==0, "xmpio_write_frmarkers") 791 #endif 792 end if 793 794 call MPI_FILE_CLOSE(unt,mpierr) 795 ABI_CHECK_MPI(mpierr,"FILE_CLOSE!") 796 797 ! Add full pawrhoij datastructure at the end of the file. 798 !if (present(pawrhoij_all) .and. me_fft == master .and. hdr%usepaw == 1) then 799 ! if (open_file(path, msg, newunit=unt, form='unformatted', status='old', action="write", access="append") /= 0) then 800 ! ABI_ERROR(msg) 801 ! end if 802 ! call pawrhoij_io(pawrhoij_all,un,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Write") 803 ! close(unt) 804 !end if 805 #endif 806 807 case (IO_MODE_ETSF) 808 file_etsf = nctk_ncify(path) 809 810 ! Write datar. 811 ncerr = nctk_write_datar(varname,file_etsf,ngfft,cplex,nfft,nspden, & 812 comm_fft,fftn3_distrib,ffti3_local,datar,action="create") 813 NCF_CHECK(ncerr) 814 call xmpi_barrier(comm_fft) 815 816 ! Master writes the header. 817 if (xmpi_comm_rank(comm_fft) == master) then 818 NCF_CHECK(nctk_open_modify(ncid, file_etsf, xmpi_comm_self)) 819 NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.)) 820 ! Add information on the crystalline structure. 821 NCF_CHECK(crystal%ncwrite(ncid)) 822 if (present(ebands)) then 823 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 824 end if 825 826 ! Add full pawrhoij datastructure. 827 !if (present(pawrhoij_all) .and. me_fft == master .and. hdr%usepaw == 1) then 828 ! call pawrhoij_io(pawrhoij_all,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Write", form="netcdf") 829 !end if 830 831 NCF_CHECK(nf90_close(ncid)) 832 end if 833 834 case default 835 ABI_ERROR(sjoin("Wrong iomode:",itoa(my_iomode))) 836 end select 837 838 call cwtime_report(" IO operation", cputime, walltime, gflops) 839 840 return 841 842 ! Handle Fortran IO error 843 10 continue 844 ABI_ERROR(errmsg) 845 846 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
871 subroutine fftdatar_write_from_hdr(varname,path,iomode,hdr,ngfft,cplex,nfft,nspden,datar,mpi_enreg,eigen) 872 873 !Arguments ------------------------------------ 874 !scalars 875 integer,intent(in) :: iomode,cplex,nfft,nspden 876 character(len=*),intent(in) :: varname,path 877 type(hdr_type),intent(inout) :: hdr 878 type(MPI_type),intent(in) :: mpi_enreg 879 !arrays 880 integer,intent(in) :: ngfft(18) 881 real(dp),intent(inout) :: datar(cplex*nfft,nspden) 882 real(dp),optional,intent(in) :: eigen(:) 883 884 !Local variables------------------------------- 885 !!scalars 886 integer :: mband 887 type(crystal_t) :: crystal 888 type(ebands_t) :: ebands 889 !arrays 890 real(dp),allocatable :: ene3d(:,:,:) 891 892 ! ************************************************************************* 893 894 crystal = hdr%get_crystal() 895 896 if (present(eigen)) then 897 mband = maxval(hdr%nband) 898 ABI_CHECK(size(eigen) == mband * hdr%nkpt * hdr%nsppol, "Wrong size(eigen)") 899 ABI_MALLOC(ene3d, (mband, hdr%nkpt, hdr%nsppol)) 900 call unpack_eneocc(hdr%nkpt, hdr%nsppol, mband, hdr%nband, eigen, ene3d) 901 ebands = ebands_from_hdr(hdr, mband, ene3d) 902 ABI_FREE(ene3d) 903 904 call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg,ebands=ebands) 905 call ebands_free(ebands) 906 else 907 call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg) 908 end if 909 910 call crystal%free() 911 912 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
1253 integer function fort_denpot_skip(unit, msg) result(ierr) 1254 1255 !Arguments ------------------------------------ 1256 integer,intent(in) :: unit 1257 character(len=*),intent(out) :: msg 1258 1259 !Local variables------------------------------- 1260 integer :: ii,fform,nspden 1261 type(hdr_type) :: hdr 1262 1263 ! ********************************************************************* 1264 1265 ierr = 1 1266 call hdr_fort_read(hdr, unit, fform) 1267 if (fform == 0) then 1268 msg = "hdr_fort_read returned fform == 0"; return 1269 end if 1270 1271 nspden = hdr%nspden 1272 call hdr%free() 1273 1274 ! Skip the records with v1. 1275 do ii=1,nspden 1276 read(unit, iostat=ierr, iomsg=msg) 1277 if (ierr /= 0) return 1278 end do 1279 1280 ierr = 0 1281 1282 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),target :: 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 type(pawrhoij_type),pointer:: pawrhoij__(:) 175 176 ! ************************************************************************* 177 178 DBG_ENTER("COLL") 179 180 ncplxfft = cplex*nfft 181 182 restartpaw=0 183 my_fildata = fildata 184 nspden = dtset%nspden; usewvl = dtset%usewvl 185 186 ! Check validity of arguments--only rho(r) (51,52) and V(r) (101,102) are presently supported 187 if ( (fform-1)/2 /=25 .and. (fform-1)/2 /=50 ) then 188 write(msg,'(a,i0,a)')' Input fform= ',fform,' not allowed.' 189 ABI_BUG(msg) 190 end if 191 192 ! Print input fform 193 if ( (fform-1)/2==25 .and. rdwr==1) then 194 msg = ' ioarr: reading density data ' 195 else if ( (fform-1)/2==25 .and. rdwr==2) then 196 msg = ' ioarr: writing density data' 197 else if ( (fform-1)/2==50 .and. rdwr==1) then 198 msg = ' ioarr: reading potential data' 199 else if ( (fform-1)/2==50 .and. rdwr==2) then 200 msg = ' ioarr: writing potential data' 201 end if 202 call wrtout(std_out,msg) 203 204 call wrtout(std_out, 'ioarr: file name is: '//TRIM(fildata)) 205 206 if (accessfil == IO_MODE_ETSF) then ! Initialize filename in case of ETSF file. 207 file_etsf = nctk_ncify(fildata) 208 call wrtout(std_out,sjoin('file name for ETSF access: ', file_etsf)) 209 end if 210 211 !Some definitions for MPI-IO access 212 spaceComm = mpi_enreg%comm_cell 213 comm_cell = mpi_enreg%comm_cell 214 comm_fft = mpi_enreg%comm_fft 215 216 if (accessfil == 4) then 217 iomode=IO_MODE_MPI 218 if (rdwr==1) then 219 spaceComm=mpi_enreg%comm_cell 220 else 221 spaceComm=mpi_enreg%comm_fft 222 end if 223 me=xmpi_comm_rank(spaceComm) 224 if (mpi_enreg%nproc_fft>1) then 225 me_fft=mpi_enreg%me_fft 226 spaceComm_io=mpi_enreg%comm_fft 227 else 228 me_fft=0 229 spaceComm_io=xmpi_comm_self 230 end if 231 end if 232 if (usewvl==1) then 233 spaceComm=mpi_enreg%comm_cell 234 me=xmpi_comm_rank(spaceComm) 235 end if 236 237 ! Change communicators and ranks if we are calling ioarr with one single processor. 238 if (present(single_proc)) then 239 if (single_proc) then 240 spaceComm = xmpi_comm_self 241 spaceComm_io = xmpi_comm_self 242 ABI_CHECK(mpi_enreg%nproc_fft == 1, "single_proc cannot be used when nproc_fft > 1") 243 comm_cell = xmpi_comm_self 244 comm_fft = xmpi_comm_self 245 me = 0 246 end if 247 end if 248 249 !======================================= 250 !Handle input from disk file 251 !======================================= 252 253 call cwtime(cputime, walltime, gflops, "start") 254 255 if (rdwr==1) then 256 if (accessfil == 0 .or. accessfil == 4) then 257 258 ! Here master checks if the input rho(r) is given on a FFT mesh that quals 259 ! the one used in the run. If not, we perform a Fourier interpolation, we write the 260 ! interpolated rho(r) to a temporary file and we use this file to restart. 261 if (ALLOW_FFTINTERP .and. usewvl==0) then 262 need_fftinterp = .False.; icheck_fft = .True. 263 ! only master checks the FFT mesh if MPI-IO. All processors read ngfft if Fortran-IO 264 ! Note that, when Fortran-IO is used, we don't know if the routine is called 265 ! by a single processor or by all procs in comm_cell hence we cannot broadcast my_fildata 266 ! inside spaceComm as done if accessfil == 4 267 if (accessfil == 4) icheck_fft = (xmpi_comm_rank(spaceComm)==master) 268 269 if (icheck_fft) then 270 if (open_file(fildata,msg,newunit=in_unt,form='unformatted',status='old') /= 0) then 271 ABI_ERROR(msg) 272 end if 273 274 call hdr_io(fform_dum,hdr0,rdwr,in_unt) 275 need_fftinterp = (ANY(hdr%ngfft/=hdr0%ngfft) ) 276 qeq0=(hdr%qptn(1)**2+hdr%qptn(2)**2+hdr%qptn(3)**2<1.d-14) 277 ! FIXME: SHould handle double-grid if PAW 278 nfftot_in = product(hdr0%ngfft(1:3)) 279 nfftot_out = product(hdr%ngfft(1:3)) 280 281 if (need_fftinterp) then 282 write(msg, "(2a,2(a,3(i0,1x)))")& 283 "Will perform Fourier interpolation since in and out ngfft differ",ch10,& 284 "ngfft in file: ",hdr0%ngfft,", expected ngfft: ",hdr%ngfft 285 ABI_WARNING(msg) 286 287 ! Read rho(r) from file, interpolate it, write data and change fildata 288 ABI_MALLOC(rhor_file, (cplex*nfftot_in, hdr0%nspden)) 289 ABI_MALLOC(rhog_in, (2, nfftot_in)) 290 ABI_MALLOC(rhor_out, (cplex*nfftot_out, hdr0%nspden)) 291 ABI_MALLOC(rhog_out, (2, nfftot_out)) 292 293 do ispden=1,hdr0%nspden 294 read(in_unt, err=10, iomsg=errmsg) (rhor_file(iarr,ispden), iarr=1,cplex*nfftot_in) 295 end do 296 297 ngfft_in = dtset%ngfft; ngfft_out = dtset%ngfft 298 ngfft_in(1:3) = hdr0%ngfft(1:3); ngfft_out(1:3) = hdr%ngfft(1:3) 299 ngfft_in(4:6) = hdr0%ngfft(1:3); ngfft_out(4:6) = hdr%ngfft(1:3) 300 ngfft_in(9:18) = 0; ngfft_out(9:18) = 0 301 ngfft_in(10) = 1; ngfft_out(10) = 1 302 303 call initmpi_seq(MPI_enreg_seq) 304 ! Which one is coarse? Note that this part is not very robust and can fail! 305 if (ngfft_in(2) * ngfft_in(3) < ngfft_out(2) * ngfft_out(3)) then 306 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft_in(2),ngfft_in(3),'all') 307 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfft_out(2),ngfft_out(3),'all') 308 else 309 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfft_in(2),ngfft_in(3),'all') 310 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft_out(2),ngfft_out(3),'all') 311 end if 312 313 call fourier_interpol(cplex,hdr0%nspden,0,0,nfftot_in,ngfft_in,nfftot_out,ngfft_out,& 314 MPI_enreg_seq,rhor_file,rhor_out,rhog_in,rhog_out) 315 316 call destroy_mpi_enreg(MPI_enreg_seq) 317 318 ! MG Hack: Change fildata so that we will use this file to read the correct rho(r) 319 ! FIXME: This should be done in a cleaner way! 320 my_fildata = trim(fildata)//"__fftinterp_rhor__" 321 if (my_fildata == fildata) my_fildata = "__fftinterp_rhor__" 322 if (open_file(my_fildata,msg,newunit=out_unt,form='unformatted',status='unknown') /= 0) then 323 ABI_ERROR(msg) 324 end if 325 call hdr_io(fform_dum,hdr,2,out_unt) 326 do ispden=1,hdr0%nspden 327 write(out_unt, err=10, iomsg=errmsg) (rhor_out(iarr,ispden),iarr=1,cplex*nfftot_out) 328 end do 329 close(out_unt) 330 331 ABI_FREE(rhor_file) 332 ABI_FREE(rhog_in) 333 ABI_FREE(rhor_out) 334 ABI_FREE(rhog_out) 335 end if ! need_fftinterp 336 337 call hdr0%free() 338 close(in_unt, err=10, iomsg=errmsg) 339 end if ! master 340 if (accessfil == 4) call xmpi_bcast(my_fildata,master,spaceComm,ierr) 341 end if 342 343 if (accessfil == 4) then 344 unt = get_unit() 345 call WffOpen(iomode,spaceComm,my_fildata,ierr,wff,0,me,unt,spaceComm_io) 346 call hdr_io(fform_dum,hdr0,rdwr,wff) 347 ! Compare the internal header and the header from the file 348 call hdr_check(fform,fform_dum,hdr,hdr0,'COLL',restart,restartpaw) 349 350 else 351 if (open_file(my_fildata, msg, newunit=unt, form="unformatted", status="old", action="read") /= 0) then 352 ABI_ERROR(msg) 353 end if 354 ! Initialize hdr0, thanks to reading of unwff1 355 call hdr_io(fform_dum,hdr0,rdwr,unt) 356 ! Compare the internal header and the header from the file 357 call hdr_check(fform,fform_dum,hdr,hdr0,'COLL',restart,restartpaw) 358 end if 359 etotal=hdr0%etot 360 361 ! NOTE: should check that restart is possible !! 362 !call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 363 364 ! If nspden[file] /= nspden, need a temporary array 365 if (hdr0%nspden/=nspden) then 366 ABI_MALLOC(arr_file,(cplex*nfft,hdr0%nspden)) 367 else 368 arr_file => arr 369 end if 370 371 ! Read data 372 do ispden=1,hdr0%nspden 373 if(accessfil == 4) then 374 call xderiveRRecInit(wff,ierr) 375 call xderiveRead(wff,arr(1:ncplxfft,ispden),ncplxfft,spaceComm_io,ierr) 376 call xderiveRRecEnd(wff,ierr) 377 378 !call xmpio_read_dp(mpi_fh,offset,sc_mode,ncount,buf,fmarker,mpierr,advance) 379 !do idat=1,ndat 380 ! do i3=1,n3 381 ! if( fftn3_distrib(i3) == me_fft) then 382 ! i3_local = ffti3_local(i3) 383 ! i3_ldat = i3_local + (idat - 1) * nd3proc 384 ! do i2=1,n2 385 ! frbase=n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * nfft 386 ! do i1=1,n1 387 ! fofr(i1+frbase)=workr(1,i1,i2,i3_ldat) 388 ! end do 389 ! end do 390 ! end if 391 ! end do 392 !end do 393 394 else 395 read(unt, err=10, iomsg=errmsg) (arr_file(iarr,ispden),iarr=1,ncplxfft) 396 end if 397 end do 398 399 if (accessfil == 4) then 400 call wffclose(wff,ierr) 401 else 402 close (unit=unt, err=10, iomsg=errmsg) 403 end if 404 405 else if (accessfil == 3) then 406 407 ! Read the header and broadcast it in comm_cell 408 ! FIXME: Use xmpi_comm_self for the time-being because, in loper, ioarr 409 ! is called by me==0 410 call hdr_read_from_fname(hdr0, file_etsf, fform_dum, comm_cell) 411 !call hdr_read_from_fname(hdr0, file_etsf, fform_dum, xmpi_comm_self) 412 ABI_CHECK(fform_dum/=0, "hdr_read_from_fname returned fform 0") 413 414 ! Compare the internal header and the header from the file 415 call hdr_check(fform, fform_dum, hdr, hdr0, 'COLL', restart, restartpaw) 416 417 ! If nspden[file] /= nspden, need a temporary array 418 if (hdr0%nspden /= nspden) then 419 ABI_MALLOC(arr_file,(cplex*nfft,hdr0%nspden)) 420 else 421 arr_file => arr 422 end if 423 424 if (usewvl == 1) then 425 ! Read the array 426 if (fform==52) then ! density 427 varname = "density" 428 else if (fform==102) then ! all potential forms!!!! 429 varname = "exchange_correlation_potential" 430 end if 431 432 ! Open the file 433 NCF_CHECK(nctk_open_read(ncid, file_etsf, xmpi_comm_self)) 434 NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, varname), arr_file)) 435 NCF_CHECK(nf90_close(ncid)) 436 else 437 ! Get MPI-FFT tables from input ngfft 438 call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 439 440 ! Get the name of the netcdf variable from the ABINIT extension and read data. 441 varname = varname_from_fname(file_etsf) 442 ncerr = nctk_read_datar(file_etsf,varname,ngfft,cplex,nfft,hdr0%nspden,comm_fft,fftn3_distrib,ffti3_local,arr) 443 NCF_CHECK(ncerr) 444 end if 445 446 else 447 write(msg,'(a,i0,a)')'Bad value for accessfil', accessfil, ' on read ' 448 ABI_BUG(msg) 449 end if 450 451 call wrtout(std_out,sjoin("data read from disk file: ", fildata)) 452 453 etotal=hdr0%etot 454 455 ! Possibly need to convert the potential/density spin components 456 if (hdr0%nspden/=nspden) then 457 call denpot_spin_convert(arr_file,hdr0%nspden,arr,nspden,fform) 458 ABI_FREE(arr_file) 459 end if 460 461 ! Eventually copy (or distribute) PAW data 462 if (rdwrpaw==1.and.restartpaw/=0) then 463 pawrhoij__ => hdr0%pawrhoij ! Trick needed by nvhpc 23.9 464 if (size(pawrhoij__) /= size(pawrhoij)) then 465 call pawrhoij_copy(hdr0%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, & 466 keep_nspden=.true.) 467 else 468 call pawrhoij_copy(hdr0%pawrhoij,pawrhoij,keep_nspden=.true.) 469 end if 470 end if 471 472 if (accessfil == 0 .or. accessfil == 3 .or. accessfil == 4) call hdr0%free() 473 474 ! ======================================= 475 ! Set up for writing data 476 ! ======================================= 477 else if (rdwr==2) then 478 479 ! In the wavelet case (isolated boundary counditions), the 480 ! arr array has a buffer that we need to remove. 481 if (usewvl == 1) then 482 #ifdef HAVE_BIGDFT 483 zindex = wvl_den%denspot%dpbox%nscatterarr(me, 3) 484 if (wvl_den%denspot%rhod%geocode == 'F') then 485 n1 = (wvl_den%denspot%dpbox%ndims(1) - 31) / 2 486 n2 = (wvl_den%denspot%dpbox%ndims(2) - 31) / 2 487 n3 = (wvl_den%denspot%dpbox%ndims(3) - 31) / 2 488 zstart = max(15 - zindex, 0) 489 zstop = wvl_den%denspot%dpbox%nscatterarr(me, 2) + & 490 & wvl_den%denspot%dpbox%nscatterarr(me, 4) - & 491 & max(zindex + wvl_den%denspot%dpbox%nscatterarr(me, 2) & 492 & - 2 * n3 - 15, 0) 493 else 494 ABI_ERROR('ioarr: WVL not implemented yet.') 495 end if 496 if (zstop - zstart + 1 > 0) then 497 ! Our slab contains (zstop - zstart + 1) elements 498 ABI_MALLOC(my_density,((n1*2)*(n2*2)*(zstop-zstart),nspden)) 499 ! We copy the data except the buffer to my_density 500 ind = 0 501 502 do i3 = zstart, zstop - 1, 1 503 ia = (i3 - 1) * dtset%ngfft(1) * dtset%ngfft(2) 504 do i2 = 0, 2 * n2 - 1, 1 505 i = ia + (i2 + 14) * dtset%ngfft(1) + 14 506 do i1 = 0, 2 * n1 - 1, 1 507 i = i + 1 508 ind = ind + 1 509 my_density(ind, :) = arr(i, :) 510 end do 511 end do 512 end do 513 else 514 nullify(my_density) 515 end if 516 #else 517 BIGDFT_NOTENABLED_ERROR() 518 if(.false. .and. present(wvl_den))then 519 write(std_out,*)' One should not be here' 520 endif 521 #endif 522 end if 523 524 ! Make sure ngfft agrees with hdr%ngfft. 525 if (usewvl == 0) then 526 if (any(ngfft(:3) /= hdr%ngfft(:3))) then 527 write(msg,"(2(a,3(1x,i0)))")"input ngfft: ",ngfft(:3),"differs from hdr%ngfft: ",hdr%ngfft(:3) 528 ABI_ERROR(msg) 529 end if 530 end if 531 532 if (accessfil == 0 .or. accessfil == 4) then 533 if(accessfil == 4) then 534 unt = get_unit() 535 call WffOpen(iomode,spaceComm,fildata,ierr,wff,0,me,unt) 536 call hdr_io(fform,hdr,rdwr,wff) 537 else 538 if (open_file(fildata, msg, newunit=unt, form='unformatted', status='unknown', action="write") /= 0) then 539 ABI_ERROR(msg) 540 end if 541 542 ! Write header 543 call hdr_io(fform,hdr,rdwr,unt) 544 end if 545 546 ! Write actual data 547 do ispden=1,nspden 548 if(accessfil == 4) then 549 call xderiveWRecInit(wff,ierr,me_fft) 550 call xderiveWrite(wff,arr(1:ncplxfft,ispden),ncplxfft,spaceComm_io,ierr) 551 call xderiveWRecEnd(wff,ierr,me_fft) 552 else 553 if (usewvl == 0) then 554 write(unt, err=10, iomsg=errmsg) (arr(iarr,ispden),iarr=1,ncplxfft) 555 else 556 write(unt, err=10, iomsg=errmsg) (my_density(iarr,ispden),iarr=1,size(my_density, 1)) 557 end if 558 end if 559 end do 560 561 if(accessfil == 4) then 562 call WffClose(wff,ierr) 563 else 564 close(unt, err=10, iomsg=errmsg) 565 end if 566 567 else if ( accessfil == 3 ) then 568 569 ! Master in comm_fft creates the file and writes the header. 570 if (xmpi_comm_rank(comm_fft) == 0) then 571 call hdr%write_to_fname(file_etsf, fform) 572 end if 573 call xmpi_barrier(comm_fft) 574 575 ! Write the array 576 if (usewvl == 0) then 577 ! Get MPI-FFT tables from input ngfft 578 call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 579 580 varname = varname_from_fname(file_etsf) 581 ncerr = nctk_write_datar(varname,file_etsf,ngfft,cplex,nfft,nspden,comm_fft,fftn3_distrib,ffti3_local,arr) 582 NCF_CHECK(ncerr) 583 else 584 NCF_CHECK(nctk_open_modify(ncid, file_etsf, xmpi_comm_self)) 585 586 if (fform==52) then ! density 587 varname = "density" 588 if (usewvl == 0) then 589 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), arr)) 590 else 591 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), my_density)) 592 end if 593 else if (fform==102) then ! all potential forms!!!! 594 varname = "exchange_correlation_potential" 595 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), arr)) 596 end if 597 598 NCF_CHECK(nf90_close(ncid)) 599 end if 600 601 else 602 write(msg,'(a,i0,a)')'Bad value for accessfil', accessfil, ' on write ' 603 ABI_ERROR(msg) 604 end if 605 606 if (usewvl == 1 .and. associated(my_density)) then 607 ABI_FREE(my_density) 608 end if 609 610 call wrtout(std_out,sjoin(' Data written to disk file:', fildata)) 611 612 else 613 write(msg,'(a,i0,a)')'Called with rdwr = ',rdwr,' not allowed.' 614 ABI_BUG(msg) 615 end if 616 617 call cwtime_report(" IO operation", cputime, walltime, gflops) 618 619 DBG_EXIT("COLL") 620 621 return 622 623 ! Handle Fortran IO error 624 10 continue 625 ABI_ERROR(errmsg) 626 627 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
965 subroutine read_rhor(fname, cplex, nspden, nfft, ngfft, pawread, mpi_enreg, orhor, ohdr, pawrhoij, comm, & 966 check_hdr, allow_interp) ! Optional 967 968 !Arguments ------------------------------------ 969 !scalars 970 integer,intent(in) :: cplex,nfft,nspden,pawread,comm 971 character(len=*),intent(in) :: fname 972 type(MPI_type),intent(in) :: mpi_enreg 973 type(hdr_type),intent(out) :: ohdr 974 type(hdr_type),optional,intent(in) :: check_hdr 975 logical,optional,intent(in) :: allow_interp 976 !arrays 977 integer,intent(in) :: ngfft(18) 978 real(dp),intent(out) :: orhor(cplex*nfft,nspden) 979 type(pawrhoij_type),intent(inout) :: pawrhoij(:) 980 981 !Local variables------------------------------- 982 !scalars 983 integer,parameter :: master=0 984 integer :: unt,fform,iomode,my_rank,mybase,globase,cplex_file 985 integer :: ispden,ifft,nfftot_file,nprocs,ierr,i1,i2,i3,i3_local,n1,n2,n3 986 integer,parameter :: fform_den=52 987 integer :: restart, restartpaw 988 integer :: ncerr 989 real(dp) :: ratio,ucvol 990 real(dp) :: cputime,walltime,gflops 991 logical :: need_interp,have_mpifft,allow_interp__ 992 character(len=500) :: msg,errmsg 993 character(len=fnlen) :: my_fname 994 character(len=nctk_slen) :: varname 995 !arrays 996 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:) 997 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2) 998 real(dp),allocatable :: rhor_file(:,:),rhor_tmp(:,:) 999 type(pawrhoij_type),allocatable :: pawrhoij_file(:) 1000 1001 ! ************************************************************************* 1002 1003 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1004 n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); have_mpifft = (nfft /= product(ngfft(1:3))) 1005 allow_interp__ = .False.; if (present(allow_interp)) allow_interp__ = allow_interp 1006 1007 call timab(1280,1,tsec) 1008 call wrtout(std_out, sjoin(" About to read data(r) from:", fname), do_flush=.True.) 1009 call cwtime(cputime, walltime, gflops, "start") 1010 1011 ! Master node opens the file, read the header and the FFT data 1012 ! This approach facilitates the interpolation of the density if in_ngfft(1:3) /= file_ngfft(1:3) 1013 if (my_rank == master) then 1014 my_fname = fname 1015 if (nctk_try_fort_or_ncfile(my_fname, msg) /= 0 ) then 1016 ABI_ERROR(msg) 1017 end if 1018 1019 iomode = iomode_from_fname(my_fname) 1020 select case (iomode) 1021 1022 case (IO_MODE_FORTRAN, IO_MODE_MPI) 1023 if (open_file(my_fname, msg, newunit=unt, form='unformatted', status='old', action="read") /= 0) then 1024 ABI_ERROR(msg) 1025 end if 1026 1027 call hdr_fort_read(ohdr, unt, fform) 1028 1029 ! Check important dimensions. 1030 ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname)) 1031 !if (fform /= fform_den) then 1032 ! write(msg, "(3a, 2(a, i0))")' File: ',trim(my_fname),ch10,' is not a density file. fform: ',fform,", expecting: ", fform_den 1033 ! ABI_WARNING(msg) 1034 !end if 1035 cplex_file = 1 1036 if (ohdr%pertcase /= 0) then 1037 cplex_file = 2; if (ohdr%qptn(1)**2 + ohdr%qptn(2)**2 + ohdr%qptn(3)**2 <1.d-14) cplex_file= 1 1038 end if 1039 ABI_CHECK(cplex_file == cplex, "cplex_file != cplex") 1040 1041 ! Read FFT array (full box) 1042 nfftot_file = product(ohdr%ngfft(:3)) 1043 ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden)) 1044 do ispden=1,ohdr%nspden 1045 read(unt, err=10, iomsg=errmsg) (rhor_file(ifft,ispden), ifft=1,cplex*nfftot_file) 1046 end do 1047 close(unt) 1048 1049 case (IO_MODE_ETSF) 1050 NCF_CHECK(nctk_open_read(unt, my_fname, xmpi_comm_self)) 1051 call hdr_ncread(ohdr, unt, fform) 1052 1053 ! Check important dimensions. 1054 ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname)) 1055 !if (fform /= fform_den) then 1056 ! write(msg, "(2a, 2(a, i0))")' File: ',trim(my_fname),' is not a density file: fform= ',fform,", expecting:", fform_den 1057 ! ABI_WARNING(msg) 1058 !end if 1059 1060 cplex_file = 1 1061 if (ohdr%pertcase /= 0) then 1062 cplex_file = 2; if (ohdr%qptn(1)**2 + ohdr%qptn(2)**2 + ohdr%qptn(3)**2 <1.d-14) cplex_file= 1 1063 end if 1064 ABI_CHECK(cplex_file == cplex, "cplex_file != cplex") 1065 1066 ! Read FFT array (full box) 1067 nfftot_file = product(ohdr%ngfft(:3)) 1068 ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden)) 1069 1070 varname = varname_from_fname(my_fname) 1071 ncerr= nf90_get_var(unt, nctk_idname(unt, varname), rhor_file, & 1072 count=[cplex, ohdr%ngfft(1), ohdr%ngfft(2), ohdr%ngfft(3), ohdr%nspden]) 1073 NCF_CHECK(ncerr) 1074 NCF_CHECK(nf90_close(unt)) 1075 1076 case default 1077 ABI_ERROR(sjoin("Wrong iomode:", itoa(iomode))) 1078 end select 1079 1080 need_interp = any(ohdr%ngfft(1:3) /= ngfft(1:3)) 1081 if (need_interp) then 1082 msg = sjoin("Different FFT meshes. Caller expects:", ltoa(ngfft(1:3)), & 1083 ". File: ", ltoa(ohdr%ngfft(1:3)), ". Need to perform interpolation.") 1084 ABI_COMMENT(msg) 1085 if (.not. allow_interp__) then 1086 write(msg, "(5a)") & 1087 " Cannot continue as allow_interp = .False. ", ch10, & 1088 " Please set ngfft to: ", trim(ltoa(ohdr%ngfft(1:3))), " in the input file" 1089 ABI_ERROR(msg) 1090 end if 1091 1092 ABI_MALLOC(rhor_tmp, (cplex*product(ngfft(1:3)), ohdr%nspden)) 1093 call timab(1281,1,tsec) 1094 call interpolate_denpot(cplex, ohdr%ngfft(1:3), ohdr%nspden, rhor_file, ngfft(1:3), rhor_tmp) 1095 call timab(1281,2,tsec) 1096 1097 ohdr%ngfft(1:3) = ngfft(1:3) 1098 nfftot_file = product(ohdr%ngfft(:3)) 1099 ABI_FREE(rhor_file) 1100 ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden)) 1101 rhor_file = rhor_tmp 1102 ABI_FREE(rhor_tmp) 1103 1104 ! Renormalize charge to avoid errors due to the interpolation. 1105 ! Do this only for NC since for PAW we should add the onsite contribution. 1106 ! This is left to the caller. 1107 !if (ohdr%usepaw == 0) then 1108 if (ohdr%usepaw == 0 .and. fform == fform_den) then 1109 call metric(gmet, gprimd, -1, rmet, ohdr%rprimd, ucvol) 1110 ratio = ohdr%nelect / (sum(rhor_file(:,1))*ucvol/ product(ngfft(1:3))) 1111 rhor_file = rhor_file * ratio 1112 write(msg,'(a,f8.2,a,f8.4)')' Expected nelect: ',ohdr%nelect,' renormalization ratio: ',ratio 1113 call wrtout(std_out,msg) 1114 end if 1115 end if ! need_interp 1116 1117 ! Read PAW Rhoij 1118 if (ohdr%usepaw == 1) then 1119 ABI_MALLOC(pawrhoij_file, (ohdr%natom)) 1120 call pawrhoij_nullify(pawrhoij_file) 1121 call pawrhoij_alloc(pawrhoij_file, ohdr%pawrhoij(1)%cplex_rhoij, ohdr%pawrhoij(1)%nspden, ohdr%pawrhoij(1)%nspinor, & 1122 ohdr%pawrhoij(1)%nsppol, ohdr%typat, lmnsize=ohdr%lmn_size, qphase=ohdr%pawrhoij(1)%qphase) 1123 call pawrhoij_copy(ohdr%pawrhoij, pawrhoij_file) 1124 end if 1125 1126 ! Check that restart is possible ! 1127 ! This check must be done here because we may have changed hdr% if need_interp 1128 if (present(check_hdr)) then 1129 ! FIXME: Temporary hack: fform_den to make hdr_check happy! 1130 call hdr_check(fform_den, fform_den, check_hdr, ohdr, "COLL", restart, restartpaw) 1131 !call hdr_check(fform_den, fform, check_hdr, ohdr, "COLL", restart, restartpaw) 1132 end if 1133 1134 end if ! master 1135 1136 if (nprocs == 1) then 1137 if (ohdr%nspden == nspden) then 1138 orhor = rhor_file 1139 else 1140 call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform) 1141 end if 1142 if (pawread == 1) call pawrhoij_copy(pawrhoij_file, pawrhoij, keep_nspden=.true.) 1143 1144 else 1145 call ohdr%bcast(master, my_rank, comm) 1146 call xmpi_bcast(fform, master, comm, ierr) 1147 1148 ! Eventually copy (or distribute) PAW data 1149 if (ohdr%usepaw == 1 .and. pawread == 1) then 1150 if (my_rank /= master) then 1151 ABI_MALLOC(pawrhoij_file, (ohdr%natom)) 1152 call pawrhoij_nullify(pawrhoij_file) 1153 call pawrhoij_alloc(pawrhoij_file, ohdr%pawrhoij(1)%cplex_rhoij, ohdr%pawrhoij(1)%nspden, ohdr%pawrhoij(1)%nspinor, & 1154 ohdr%pawrhoij(1)%nsppol, ohdr%typat, lmnsize=ohdr%lmn_size, qphase=ohdr%pawrhoij(1)%qphase) 1155 end if 1156 if (size(ohdr%pawrhoij) /= size(pawrhoij)) then 1157 call pawrhoij_copy(ohdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, & 1158 keep_nspden=.true.) 1159 else 1160 call pawrhoij_copy(ohdr%pawrhoij,pawrhoij, keep_nspden=.true.) 1161 end if 1162 end if 1163 1164 if (my_rank /= master) then 1165 nfftot_file = product(ohdr%ngfft(1:3)) 1166 ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden)) 1167 end if 1168 call xmpi_bcast(rhor_file, master, comm,ierr) 1169 1170 if (have_mpifft) then 1171 ! Extract slice treated by this MPI-FFT process. 1172 call ptabs_fourdp(mpi_enreg, ngfft(2), ngfft(3), fftn2_distrib, ffti2_local, fftn3_distrib, ffti3_local) 1173 if (ohdr%nspden==nspden) then 1174 do ispden=1,nspden 1175 do i3=1,n3 1176 if (fftn3_distrib(i3) /= mpi_enreg%me_fft) cycle 1177 i3_local = ffti3_local(i3) 1178 do i2=1,n2 1179 mybase = cplex * (n1 * (i2-1 + n2*(i3_local-1))) 1180 globase = cplex * (n1 * (i2-1 + n2*(i3-1))) 1181 do i1=1,n1*cplex 1182 orhor(i1+mybase,ispden) = rhor_file(i1+globase,ispden) 1183 end do 1184 end do 1185 end do 1186 end do 1187 else 1188 do i3=1,n3 1189 if (fftn3_distrib(i3) /= mpi_enreg%me_fft) cycle 1190 i3_local = ffti3_local(i3) 1191 do i2=1,n2 1192 mybase = 1 + cplex * (n1 * (i2-1 + n2*(i3_local-1))) 1193 globase = 1 + cplex * (n1 * (i2-1 + n2*(i3-1))) 1194 call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform,& 1195 istart_in=globase,istart_out=mybase,nelem=n1*cplex) 1196 end do 1197 end do 1198 end if 1199 else 1200 if (ohdr%nspden==nspden) then 1201 orhor = rhor_file 1202 else 1203 call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform) 1204 end if 1205 end if 1206 end if ! nprocs > 1 1207 1208 ABI_FREE(rhor_file) 1209 1210 if (allocated(pawrhoij_file)) then 1211 call pawrhoij_free(pawrhoij_file) 1212 ABI_FREE(pawrhoij_file) 1213 end if 1214 1215 ! Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities 1216 ! Add a small real to the magnetization 1217 if (nspden==4) orhor(:,4)=orhor(:,4)+tol14 1218 if (ohdr%usepaw==1.and.size(pawrhoij)>0) then 1219 if (pawrhoij(1)%nspden==4) then 1220 do i1=1,size(pawrhoij) 1221 pawrhoij(i1)%rhoijp(:,4)=pawrhoij(i1)%rhoijp(:,4)+tol10 1222 end do 1223 end if 1224 end if 1225 1226 call timab(1280,2,tsec) 1227 call cwtime_report(" read_rhor", cputime, walltime, gflops) 1228 return 1229 1230 ! Handle Fortran IO error 1231 10 continue 1232 ABI_ERROR(errmsg) 1233 1234 end subroutine read_rhor