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