TABLE OF CONTENTS
- ABINIT/m_iowf
- m_iowf/cg2seqblocks
- m_iowf/cg_ncwrite
- m_iowf/kg2seqblocks
- m_iowf/ncwrite_eigen1_occ
- m_iowf/outresid
- m_iowf/outwf
- m_iowf/prtkbff
ABINIT/m_iowf [ Modules ]
NAME
m_iowf
FUNCTION
Procedures for the IO of the WFK file.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, AR, MB, MVer, MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_iowf 23 24 use defs_basis 25 use defs_wvltypes 26 use m_abicore 27 use m_errors 28 use m_dtset 29 use m_xmpi 30 use m_wffile 31 use m_abi_etsf 32 use m_nctk 33 use m_wfk 34 #ifdef HAVE_NETCDF 35 use netcdf 36 #endif 37 use m_hdr 38 use m_ebands 39 40 use defs_abitypes, only : MPI_type 41 use m_time, only : cwtime, cwtime_report, timab 42 use m_io_tools, only : get_unit, flush_unit, iomode2str 43 use m_fstrings, only : endswith, sjoin 44 use m_numeric_tools, only : mask2blocks 45 use defs_datatypes, only : ebands_t, pseudopotential_type 46 use m_cgtools, only : cg_zcopy 47 use m_crystal, only : crystal_t 48 use m_rwwf, only : rwwf 49 use m_mpinfo, only : proc_distrb_cycle 50 use m_vkbr, only : calc_vkb 51 use m_wvl_rwwf, only : wvl_write 52 53 implicit none 54 55 private 56 57 public :: outwf 58 public :: outresid 59 public :: prtkbff ! Write KB form factors to WFK in netcdf format.
m_iowf/cg2seqblocks [ Functions ]
[ Top ] [ m_iowf ] [ Functions ]
NAME
cg2seqblocks
FUNCTION
INPUTS
npwtot_k npw_k= cg_k(2,npw_k*nband_k) gmpi2seq comm_fft=FFT communicator
OUTPUT
bstart bcount
SOURCE
1578 subroutine cg2seqblocks(npwtot_k,npw_k,nband,cg_k,gmpi2seq,comm_bandfft,bstart,bcount,my_cgblock) 1579 1580 !Arguments ------------------------------------ 1581 !scalars 1582 integer,intent(in) :: npwtot_k,npw_k,nband,comm_bandfft 1583 integer,intent(out) :: bstart,bcount 1584 !arrays 1585 integer,intent(in) :: gmpi2seq(npw_k) 1586 real(dp),intent(in) :: cg_k(2,npw_k,nband) 1587 real(dp),allocatable,intent(out) :: my_cgblock(:,:,:) 1588 1589 !Local variables------------------------------- 1590 !scalars 1591 integer :: me,nprocs,rank,ig,ib,igseq,ierr,nbb,nbmax,band 1592 !arrays 1593 integer,allocatable :: bstart_rank(:) 1594 real(dp),allocatable :: cgbuf(:,:,:) 1595 1596 ! ************************************************************************* 1597 1598 me = xmpi_comm_rank(comm_bandfft); nprocs = xmpi_comm_size(comm_bandfft) 1599 1600 ! Handle sequential case. 1601 if (nprocs == 1) then 1602 bstart = 1; bcount = nband 1603 ABI_MALLOC_OR_DIE(my_cgblock, (2, npwtot_k, nband), ierr) 1604 my_cgblock = cg_k 1605 return ! DOH 1606 end if 1607 1608 ABI_MALLOC(bstart_rank, (0:nprocs)) 1609 bstart_rank = [(1 + (nband/nprocs)*rank, rank=0,nprocs-1), 1 + nband] 1610 1611 ! Allocate MPI buffer (same size on each MPI proc) 1612 nbmax = 0 1613 do rank=0,nprocs-1 1614 nbmax = max(nbmax, bstart_rank(rank+1) - bstart_rank(rank)) 1615 end do 1616 ABI_MALLOC_OR_DIE(cgbuf, (2, npwtot_k, nbmax), ierr) 1617 1618 nbb = bstart_rank(me+1) - bstart_rank(me) 1619 ABI_MALLOC_OR_DIE(my_cgblock, (2, npwtot_k, nbb), ierr) 1620 1621 ! TODO: This should be replaced by gatherv but premature optimization.... 1622 do rank=0,nprocs-1 1623 bstart = bstart_rank(rank) 1624 bcount = bstart_rank(rank+1) - bstart_rank(rank) 1625 1626 do band=bstart, bstart+bcount-1 1627 ib = band - bstart + 1 1628 cgbuf(:,:,ib) = zero 1629 do ig=1,npw_k 1630 igseq = gmpi2seq(ig) 1631 cgbuf(:,igseq,ib) = cg_k(:,ig,band) 1632 end do 1633 end do ! band 1634 1635 call xmpi_sum_master(cgbuf,rank,comm_bandfft,ierr) 1636 if (me == rank) my_cgblock = cgbuf(:,:,:bcount) 1637 end do ! rank 1638 1639 bstart = bstart_rank(me) 1640 bcount = bstart_rank(me+1) - bstart_rank(me) 1641 1642 ABI_FREE(cgbuf) 1643 ABI_FREE(bstart_rank) 1644 1645 end subroutine cg2seqblocks
m_iowf/cg_ncwrite [ Functions ]
[ Top ] [ m_iowf ] [ Functions ]
NAME
cg_ncwrite
FUNCTION
Conduct output of a "wave-functions" file with netcdf
INPUTS
fname= character string giving the root to form the name of the output WFK or WFQ file if response==0, otherwise it is the filename. dtset <type(dataset_type)>=all input variables for this dataset hdr <type(hdr_type)>=the header of wf, den and pot files response: if == 0, GS wavefunctions , if == 1, RF wavefunctions mpw=maximum number of plane waves mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mkmem=maximum number of k-points treated by this node nkpt=number of k points nsppol=1 for unpolarized, 2 for spin-polarized eigen((2*mband)**response *mband*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point occ(mband*nkpt*nsppol)=occupations for all bands at each k point cg(2,mcg)=wavefunction array (storage if nkpt>1) npwarr(nkpt)=number of planewaves in basis and on boundary for each k kg(3,mpw*mkmem)=reduced planewave coordinates. mpi_enreg=information about MPI parallelization
OUTPUT
done=.True if cg_ncwrite can handle the output of the WFK file in parallel.
SOURCE
662 #ifdef HAVE_NETCDF 663 664 subroutine cg_ncwrite(fname,hdr,dtset,response,mpw,mband,nband,nkpt,nsppol,nspinor,mcg,& 665 mkmem,eigen,occ,cg,npwarr,kg,mpi_enreg,done) 666 667 !Arguments ------------------------------------ 668 !scalars 669 integer,intent(in) :: response,mband,mcg,mkmem,mpw,nkpt,nsppol,nspinor 670 character(len=*),intent(in) :: fname 671 logical,intent(out) :: done 672 type(dataset_type),intent(in) :: dtset 673 type(MPI_type),intent(in) :: mpi_enreg 674 type(hdr_type),intent(in) :: hdr 675 !arrays 676 integer, intent(in) :: nband(nkpt*nsppol),kg(3,mpw*mkmem),npwarr(nkpt) 677 real(dp),intent(in) :: cg(2,mcg) 678 real(dp),intent(in) :: eigen((2*mband)**response*mband*nkpt*nsppol),occ(mband*nkpt*nsppol) 679 680 !Local variables------------------------------- 681 !scalars 682 integer,parameter :: master=0,fform2=2 683 integer :: ii,iomode,icg,iband,ikg,ikpt,spin,me_cell,me_kpt,me_band,me_spinor,my_nspinor,nband_k,npw_k 684 integer :: comm_cell,comm_fft,comm_bandfft,mtag,formeig 685 integer :: cnt,min_cnt,max_cnt,ierr,action,source,ncid,ncerr,cg_varid,kg_varid !,eig_varid, 686 integer :: paral_kgb,npwtot_k !,start_pwblock !,start_cgblock !count_pwblock, 687 integer :: ipw,ispinor_index,npwso,npwsotot,npwtot,nspinortot,ikpt_this_proc,ispinor 688 integer :: bandpp,nproc_band,nproc_fft,nproc_spinor,me_fft,nproc_cell,nwrites 689 integer :: comm_mpiio,nranks,bstart,bcount !nbdblock,nblocks, 690 !integer :: band_blocksize,band 691 real(dp) :: cpu,wall,gflops 692 logical :: ihave_data,iam_master,single_writer,same_layout,use_collective 693 character(len=500) :: msg 694 character(len=fnlen) :: path 695 type(wfk_t) :: wfk 696 type(crystal_t) :: crystal 697 type(ebands_t) :: gs_ebands 698 !arrays 699 integer,allocatable :: kg_k(:,:),iter2kscgkg(:,:),ind_cg_mpi_to_seq(:),rank_has_cg(:),ranks_io(:)!,gblock(:,:) 700 real(dp) :: tsec(2) 701 real(dp),allocatable :: eigen3d(:,:,:),occ3d(:,:,:),cg_k(:,:),my_cgblock(:,:,:) 702 703 ! ************************************************************************* 704 705 DBG_ENTER("COLL") 706 done = .False. 707 708 path = nctk_ncify(fname) 709 call wrtout(std_out, sjoin("in cg_ncwrite with path:", path)) 710 711 ! communicators and ranks 712 comm_cell = mpi_enreg%comm_cell; me_cell = mpi_enreg%me_cell; nproc_cell = mpi_enreg%nproc_cell 713 comm_fft = mpi_enreg%comm_fft; me_fft = mpi_enreg%me_fft; nproc_fft = mpi_enreg%nproc_fft 714 comm_bandfft = mpi_enreg%comm_bandfft 715 me_kpt = mpi_enreg%me_kpt; me_band = mpi_enreg%me_band; me_spinor = mpi_enreg%me_spinor 716 iam_master = (me_kpt == master) 717 718 paral_kgb = dtset%paral_kgb 719 nproc_band = mpi_enreg%nproc_band 720 bandpp = mpi_enreg%bandpp 721 nproc_spinor = mpi_enreg%nproc_spinor 722 723 ! FIXME 724 my_nspinor = max(1, nspinor/nproc_spinor) 725 if (nspinor == 2 .and. my_nspinor == 1) then 726 ABI_ERROR("Spinor parallelization not coded yet") 727 end if 728 729 ! TODO: Be careful with response == 1 in parallel because the distribution of the cg 730 ! can be **VERY** different from cg if nprocs > nkpt * nsppol 731 !ABI_CHECK(response==0, "response == 1 not coded") 732 if (size(hdr%nband) == size(nband)) then 733 ABI_CHECK(all(Hdr%nband == nband),"nband") 734 else 735 ABI_ERROR("hdr%nband and nband have different size!") 736 end if 737 738 if (xmpi_comm_size(comm_cell) == 1) then 739 ABI_CHECK(all(npwarr == hdr%npwarr), "npwarr != hdr%npwarr") 740 end if 741 742 crystal = hdr%get_crystal() 743 744 ! TODO 745 ! Be careful with response == 1. 746 ! gs_ebands contains the GS eigenvalues and occupation and will be written if this is a 747 ! GS wfk. If we have a DFPT file, eigen stores the GKK matrix element, in this case 748 ! we don't write gs_ebands but we define new variables in the netcdf file to store the GKK 749 ABI_MALLOC(occ3d, (mband,nkpt,nsppol)) 750 call unpack_eneocc(nkpt,nsppol,mband,nband,occ,occ3d) 751 752 if (response == 0) then 753 formeig = 0 754 ABI_MALLOC(eigen3d, (mband,nkpt,nsppol)) 755 call unpack_eneocc(nkpt,nsppol,mband,nband,eigen,eigen3d) 756 gs_ebands = ebands_from_hdr(hdr, mband, eigen3d); gs_ebands%occ = occ3d 757 ABI_FREE(eigen3d) 758 else 759 formeig = 1 760 end if 761 762 ! same_layout is set to True if the internal representation of the cgs 763 ! is compatible with the representation on file. 764 iomode = IO_MODE_ETSF 765 if (response == 0) then 766 same_layout = (paral_kgb == 0 .or. (paral_kgb == 1 .and. all([nproc_fft, nproc_band, nproc_spinor] == 1))) 767 else 768 ! For the time being, these cases are not implemented in the DFPT part. 769 ABI_CHECK(nproc_fft==1, "nproc_fft != 1 not coded") 770 ABI_CHECK(nproc_band==1, "nproc_band != 1 not coded") 771 ABI_CHECK(nproc_spinor==1, "nproc_spinor != 1 not coded") 772 773 ! Note: It would be possible to use collective IO if the cg1 are block-distributed 774 same_layout = .True. 775 spin_loop: do spin=1,nsppol 776 do ikpt=1,nkpt 777 nband_k = nband(ikpt + (spin-1)*nkpt) 778 if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then 779 if (any(mpi_enreg%proc_distrb(ikpt,:nband_k,spin) /= me_kpt)) then 780 same_layout = .False.; exit spin_loop 781 end if 782 end if 783 end do 784 end do spin_loop 785 call xmpi_land(same_layout, comm_cell) 786 end if 787 788 call cwtime(cpu, wall, gflops, "start") 789 790 if (same_layout) then 791 single_writer = .True. 792 if (nctk_has_mpiio) single_writer = .False. 793 !single_writer = .True. 794 795 write(msg,'(5a,l1)')& 796 "same layout Writing WFK file: ",trim(path),", with iomode: ",trim(iomode2str(iomode)),", single writer: ",single_writer 797 call wrtout(std_out,msg,'PERS', do_flush=.True.) 798 799 if (.not. single_writer) then 800 801 #ifdef HAVE_NETCDF 802 ! master opens the file and write the metadata. 803 if (xmpi_comm_rank(comm_cell) == master) then 804 ncerr = nf90_einval 805 #ifdef HAVE_NETCDF_MPI 806 ncerr = nf90_create(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), & 807 comm=xmpi_comm_self, info=xmpio_info, ncid=ncid) 808 #endif 809 NCF_CHECK_MSG(ncerr, sjoin("create_par: ", path)) 810 811 call wfk_ncdef_dims_vars(ncid, hdr, fform2, write_hdr=.True.) 812 NCF_CHECK(crystal%ncwrite(ncid)) 813 814 if (response == 0) then 815 ! Write Gs bands 816 NCF_CHECK(ebands_ncwrite(gs_ebands, ncid)) 817 else 818 ! Write H1 matrix elements and occupancies. 819 ! Note that GS eigens are not written here. 820 call ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d) 821 end if 822 823 NCF_CHECK(nf90_close(ncid)) 824 end if 825 826 ! Compute the table for collective IO. 827 ABI_CALLOC(iter2kscgkg, (4, nkpt*nsppol)) 828 cnt = 0; icg = 0 829 do spin=1,nsppol 830 ikg = 0 831 do ikpt=1,nkpt 832 nband_k = nband(ikpt + (spin-1)*nkpt) 833 npw_k = npwarr(ikpt) 834 if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then 835 ! FIXME v2[01], v3[6], with 4 procs 836 ! v6[7] and v7[68], v7[69] fail but due to a extra line with nprocs 837 ! v7[96] with np=4 seems to be more serious but it crashes also in trunk. 838 ! v2[88] is slow likely due to outscfv 839 ABI_CHECK(all(mpi_enreg%proc_distrb(ikpt,:nband_k,spin) == me_kpt), "bands are distributed") 840 cnt = cnt + 1 841 iter2kscgkg(:,cnt) = [ikpt, spin, icg, ikg] 842 icg = icg + npw_k*my_nspinor*nband_k 843 ikg = ikg + npw_k 844 end if 845 end do 846 end do 847 if (cnt == 0) then 848 write(std_out,*)"cnt == 0 for me_cell, me, me_kpt",mpi_enreg%me_cell, mpi_enreg%me, mpi_enreg%me_kpt 849 !ABI_CHECK(cnt > 0, "This processor does not have wavefunctions!") 850 end if 851 852 call xmpi_min(cnt, min_cnt, comm_cell, ierr) 853 call xmpi_max(cnt, max_cnt, comm_cell, ierr) 854 855 ! Handle idle procs, i.e. processors that do not have wavefunctions 856 ! This happens if paral_kgb == 0 and nprocs > nkpt * nsppol (Abinit does not stop anymore!) 857 comm_mpiio = comm_cell 858 859 if (min_cnt <= 0) then 860 ABI_COMMENT("Will create subcommunicator to exclude idle processors from MPI-IO collective calls") 861 ABI_CHECK(paral_kgb == 0, "paral_kgb == 1 with idle processors should never happen") 862 863 ! Prepare the call to xmpi_subcomm that will replace comm_mpiio. 864 ABI_CALLOC(rank_has_cg, (0:nproc_cell-1)) 865 do spin=1,nsppol 866 do ikpt=1,nkpt 867 nband_k = nband(ikpt + (spin-1)*nkpt) 868 if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then 869 rank_has_cg(me_kpt) = 1 870 exit 871 end if 872 end do 873 end do 874 875 call xmpi_sum(rank_has_cg,comm_cell,ierr) 876 nranks = count(rank_has_cg == 1) 877 ABI_MALLOC(ranks_io, (nranks)) 878 cnt = 0 879 do ii=0,nproc_cell-1 880 if (rank_has_cg(ii) == 1) then 881 cnt = cnt + 1 882 ranks_io(cnt) = ii 883 end if 884 end do 885 !write(std_out,*)"nranks, ranks_io:", nranks, ranks_io 886 comm_mpiio = xmpi_subcomm(comm_cell, nranks, ranks_io) 887 if (.not. rank_has_cg(me_kpt) == 1) then 888 comm_mpiio = xmpi_comm_null 889 ABI_CHECK(rank_has_cg(me_kpt) == 0, "rank_has_cg must be 0 or 1") 890 end if 891 ABI_FREE(ranks_io) 892 ABI_FREE(rank_has_cg) 893 end if 894 895 ! Open the file in parallel inside comm_mpiio. 896 call xmpi_barrier(comm_cell) 897 if (comm_mpiio == xmpi_comm_null) goto 100 898 899 ncerr = nf90_einval 900 #ifdef HAVE_NETCDF_MPI 901 ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write),& 902 comm=comm_mpiio, info=xmpio_info, ncid=ncid) 903 #endif 904 NCF_CHECK_MSG(ncerr, sjoin("open_par: ", path)) 905 906 ! Use individual IO (default) for the G-vectors [3, mpw, nkpt] 907 NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid)) 908 spin = 1; ikg=0 909 do ikpt=1,nkpt 910 npw_k = npwarr(ikpt) 911 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,1,spin,me_kpt)) cycle 912 ncerr = nf90_put_var(ncid, kg_varid, kg(:, 1+ikg:), start=[1,1,ikpt], count=[3,npw_k,1]) 913 NCF_CHECK_MSG(ncerr, "put_kg_k") 914 ikg = ikg + npw_k 915 end do 916 917 NCF_CHECK(nf90_inq_varid(ncid, "coefficients_of_wavefunctions", cg_varid)) 918 919 use_collective = (response == 0) ! or (response == 1 .and. nproc_cell == 1) 920 921 if (use_collective) then 922 call wrtout(std_out,"Using collective IO for the CGs") 923 ! Use collective IO for the CGs 924 ncerr = nf90_einval 925 #ifdef HAVE_NETCDF_MPI 926 ncerr = nf90_var_par_access(ncid, cg_varid, nf90_collective) 927 #endif 928 NCF_CHECK(ncerr) 929 930 do cnt=1,max_cnt 931 ikpt = iter2kscgkg(1,cnt) 932 spin = iter2kscgkg(2,cnt) 933 icg = iter2kscgkg(3,cnt) 934 ikg = iter2kscgkg(4,cnt) 935 936 ! The array on file has shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 937 if (ikpt /= 0) then 938 nband_k = nband(ikpt + (spin-1)*nkpt) 939 npw_k = npwarr(ikpt) 940 941 !write(std_out,*)"About to write ikpt, spin, npw_k, icg: ",ikpt, spin, npw_k, icg 942 ncerr = nf90_put_var(ncid, cg_varid, cg(:, 1+icg:), start=[1,1,1,1,ikpt,spin], & 943 count=[2,npw_k,nspinor,nband_k,1,1]) 944 NCF_CHECK_MSG(ncerr, "writing cg") 945 else 946 ! This happens when nkpt * nsppol // nprocs != 0 947 ! Note that we are using collective MPI-IO hence all processors must call put_var 948 ! Here we re-write the ug(0) of the first (k-point, spin) treated by the node. 949 ikpt = iter2kscgkg(1,1) 950 spin = iter2kscgkg(2,1) 951 ncerr = nf90_put_var(ncid, cg_varid, cg, start=[1,1,1,1,ikpt,spin], count=[1,1,1,1,1,1]) 952 NCF_CHECK_MSG(ncerr, "re-writing cg") 953 end if 954 end do 955 956 else 957 call wrtout(std_out, "Using individual IO for the CGs") 958 ! Individual IO of the CGs (for debugging purposes) 959 icg = 0 960 do spin=1,nsppol 961 do ikpt=1,nkpt 962 nband_k = nband(ikpt + (spin-1)*nkpt) 963 npw_k = npwarr(ikpt) 964 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) cycle 965 ncerr = nf90_put_var(ncid, cg_varid, cg(:, 1+icg:), start=[1,1,1,1,ikpt,spin], count=[2,npw_k,nspinor,nband_k,1,1]) 966 NCF_CHECK_MSG(ncerr, "put_var") 967 icg = icg+npw_k*my_nspinor*nband_k 968 end do 969 end do 970 end if 971 972 NCF_CHECK(nf90_close(ncid)) 973 974 100 call xmpi_barrier(comm_cell) 975 ABI_FREE(iter2kscgkg) 976 977 done = .True. 978 call cwtime_report(" collective ncwrite", cpu, wall, gflops) 979 #endif 980 ! HAVE_NETCDF 981 else ! single_writer 982 if (nproc_cell > 1) then 983 ABI_WARNING("Slow version without MPI-IO support. Processors send data to master...") 984 else 985 call wrtout(std_out, "Using netcdf library without MPI-IO support") 986 end if 987 988 ABI_MALLOC(kg_k,(3,mpw)) 989 ABI_MALLOC_OR_DIE(cg_k,(2,mpw*my_nspinor*mband), ierr) 990 991 if (iam_master) then 992 call wfk%open_write(hdr,path,formeig,iomode,get_unit(),xmpi_comm_self,write_hdr=.True.) 993 994 NCF_CHECK(crystal%ncwrite(wfk%fh)) 995 !write(std_out,*)"after crystal_ncwrite" 996 997 ! Write eigenvalues and occupations (these arrays are not MPI-distributed) 998 if (response == 0) then 999 NCF_CHECK(ebands_ncwrite(gs_ebands, wfk%fh)) 1000 else 1001 call ncwrite_eigen1_occ(wfk%fh, nband, mband, nkpt, nsppol, eigen, occ3d) 1002 end if 1003 end if 1004 1005 icg = 0 1006 do spin=1,nsppol 1007 ikg = 0 1008 do ikpt=1,nkpt 1009 nband_k = nband(ikpt + (spin-1)*nkpt) 1010 npw_k = npwarr(ikpt) 1011 mtag = ikpt+(spin-1)*nkpt 1012 1013 call xmpi_barrier(comm_cell) 1014 1015 ! Transfer the wavefunctions and the g-vectors to the master processor 1016 source = minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,spin)) 1017 ihave_data = (source==me_kpt) 1018 1019 action=0 1020 if (iam_master .and. ihave_data) action=1 ! I am the master node, and I have the data in cg 1021 if (.not.iam_master.and.ihave_data) action=2 ! I am not the master, and I have the data => send to master 1022 if (iam_master.and..not.ihave_data) action=3 ! I am the master, and I receive the data 1023 1024 if (action==1) then ! Copy from kg and cg 1025 kg_k(:,1:npw_k) = kg(:,ikg+1:ikg+npw_k) 1026 call cg_zcopy(npw_k*my_nspinor*nband_k, cg(1,icg+1), cg_k) 1027 end if 1028 1029 ! Exchange data 1030 if (action==2.or.action==3) then 1031 call timab(48,1,tsec) 1032 if (action==2) then 1033 call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,source,kg_k,master,comm_cell,2*mtag+1,ierr) 1034 call xmpi_exch(cg(:,icg+1:icg+nband_k*npw_k*my_nspinor),2*nband_k*npw_k*my_nspinor,& 1035 & source,cg_k,master,comm_cell,2*mtag+2,ierr) 1036 else 1037 call xmpi_exch(kg_k,3*npw_k,source,kg_k,master,comm_cell,2*mtag+1,ierr) 1038 call xmpi_exch(cg_k,2*nband_k*npw_k*my_nspinor,source,cg_k,master,comm_cell,2*mtag+2,ierr) 1039 end if 1040 call timab(48,2,tsec) 1041 end if 1042 1043 ! Master writes this block of bands. 1044 if (iam_master) then 1045 if (response == 0) then 1046 call wfk%write_band_block([1,nband_k],ikpt,spin,xmpio_single,kg_k=kg_k,cg_k=cg_k,& 1047 eig_k=gs_ebands%eig(:,ikpt,spin),occ_k=gs_ebands%occ(:,ikpt,spin)) 1048 !write(std_out,*)"cg_k",cg_k(:,1:2) 1049 else 1050 call wfk%write_band_block([1,nband_k],ikpt,spin,xmpio_single,kg_k=kg_k,cg_k=cg_k) 1051 end if 1052 end if 1053 1054 if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt))) then 1055 icg = icg+npw_k*my_nspinor*nband_k 1056 ikg = ikg+npw_k 1057 end if 1058 1059 end do !ikpt 1060 end do !spin 1061 1062 ABI_FREE(kg_k) 1063 ABI_FREE(cg_k) 1064 1065 if (iam_master) call wfk%close() 1066 call xmpi_barrier(comm_cell) 1067 1068 done = .True. 1069 call cwtime_report(" individual ncwrite", cpu, wall, gflops) 1070 end if 1071 1072 else ! not same_layout 1073 1074 if (nctk_has_mpiio) then 1075 #ifdef HAVE_NETCDF 1076 call wrtout(std_out, & 1077 sjoin("scattered data. writing WFK file",trim(path),", with iomode: ",iomode2str(iomode)), 'PERS', do_flush=.True.) 1078 1079 ! master write the metadata. 1080 if (xmpi_comm_rank(comm_cell) == master) then 1081 ncerr = nf90_einval 1082 #ifdef HAVE_NETCDF_MPI 1083 ncerr = nf90_create(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), & 1084 comm=xmpi_comm_self, info=xmpio_info, ncid=ncid) 1085 #endif 1086 NCF_CHECK_MSG(ncerr, sjoin("create_par:", path)) 1087 1088 call wfk_ncdef_dims_vars(ncid, hdr, fform2, write_hdr=.True.) 1089 NCF_CHECK(crystal%ncwrite(ncid)) 1090 1091 ! Write eigenvalues and occupations (these arrays are not MPI-distributed) 1092 if (response == 0) then 1093 NCF_CHECK(ebands_ncwrite(gs_ebands, ncid)) 1094 else 1095 call ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d) 1096 end if 1097 1098 NCF_CHECK(nf90_close(ncid)) 1099 end if 1100 1101 ! Reopen the file inside comm_cell 1102 call xmpi_barrier(comm_cell) 1103 ncerr = nf90_einval 1104 #ifdef HAVE_NETCDF_MPI 1105 ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), & 1106 comm=comm_cell, info=xmpio_info, ncid=ncid) 1107 #endif 1108 NCF_CHECK_MSG(ncerr, sjoin("create_par:", path)) 1109 1110 ! Get var ids 1111 NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid)) 1112 NCF_CHECK(nf90_inq_varid(ncid, "coefficients_of_wavefunctions", cg_varid)) 1113 1114 ! Write the G-vectors 1115 ikg = 0 1116 do ikpt=1,nkpt 1117 npw_k = npwarr(ikpt) 1118 npwtot_k = hdr%npwarr(ikpt) 1119 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,1,1,me_kpt)) cycle 1120 if (me_spinor /= 0) cycle 1121 !write(std_out,*)"In G-vector loop ",ikpt,", with me_cell me_kpt me_band, me_spinor ",me_cell,me_kpt,me_band,me_spinor 1122 !write(std_out,*)" with ik, npw_k, npwtot_k: ",ikpt, npw_k, npwtot_k 1123 1124 ABI_MALLOC(ind_cg_mpi_to_seq, (npw_k)) 1125 if (allocated(mpi_enreg%my_kgtab)) then 1126 ikpt_this_proc = mpi_enreg%my_kpttab(ikpt) 1127 ind_cg_mpi_to_seq = mpi_enreg%my_kgtab(1:npw_k,ikpt_this_proc) 1128 else 1129 ABI_CHECK(nproc_fft==1, "nproc_fft !=1 and my_kgtab not allocated") 1130 ind_cg_mpi_to_seq(1:npw_k) = [(ipw, ipw=1,npw_k)] 1131 end if 1132 1133 ABI_CALLOC(kg_k, (3, npwtot_k)) 1134 do ipw=1,npw_k 1135 kg_k(:, ind_cg_mpi_to_seq(ipw)) = kg(:,ikg+ipw) 1136 end do 1137 call xmpi_sum_master(kg_k,master,comm_bandfft,ierr) 1138 if (xmpi_comm_rank(comm_bandfft) == master) then 1139 ncerr = nf90_put_var(ncid, kg_varid, kg_k, start=[1,1,ikpt], count=[3,npwtot_k,1]) 1140 NCF_CHECK_MSG(ncerr, "putting kg_k") 1141 end if 1142 ABI_FREE(kg_k) 1143 1144 ! gblock contains block-distributed G-vectors inside comm_fft 1145 !call kg2seqblocks(npwtot_k,npw_k,kg(:,ikg+1:),ind_cg_mpi_to_seq,comm_fft,start_pwblock,count_pwblock,gblock) 1146 !write(std_out,*)"gblock(:,2)",gblock(:,2) 1147 !ncerr = nf90_put_var(ncid, kg_varid, gblock, start=[1,start_pwblock,ikpt], count=[3,count_pwblock,1]) 1148 !NCF_CHECK_MSG(ncerr, "putting kg_k") 1149 !ABI_FREE(gblock) 1150 1151 ABI_FREE(ind_cg_mpi_to_seq) 1152 1153 ikg = ikg+npw_k 1154 end do 1155 1156 ! Write wavefunctions 1157 if (response == 1) then 1158 ! The cg1 is allocated with size mcg1=mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol (see dfpt_looppert) 1159 ! hence bands are not MPI distributed, this is the reason why we have to use the loop over bands 1160 ! and the check on nwrites 1161 icg = 0 1162 do spin=1,nsppol 1163 do ikpt=1,nkpt 1164 nband_k = nband(ikpt + (spin-1)*nkpt) 1165 npw_k = npwarr(ikpt) 1166 npwtot_k = hdr%npwarr(ikpt) 1167 nwrites = 0 1168 do iband=1,nband_k 1169 if (mpi_enreg%proc_distrb(ikpt,iband,spin)/=me_kpt) cycle 1170 ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 1171 nwrites = nwrites + 1 1172 ii = 1 + (iband-1)*npw_k*my_nspinor + icg 1173 ncerr = nf90_put_var(ncid, cg_varid, cg(1:,ii:), start=[1,1,1,iband,ikpt,spin], count=[2,npwtot_k,nspinor,1,1,1]) 1174 NCF_CHECK_MSG(ncerr, "put_var cg") 1175 end do ! iband 1176 if (nwrites /= 0) icg = icg + npw_k*my_nspinor*nband_k 1177 end do !ikpt 1178 end do !spin 1179 1180 else 1181 icg = 0 1182 do spin=1,nsppol 1183 do ikpt=1,nkpt 1184 nband_k = nband(ikpt + (spin-1)*nkpt) 1185 npw_k = npwarr(ikpt) 1186 npwtot_k = hdr%npwarr(ikpt) 1187 1188 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) cycle 1189 !write(std_out,*)"In u(g)-vector loop ",ikpt,", with me_cell me_kpt me_band, me_spinor ",me_cell,me_kpt,me_band,me_spinor 1190 !write(std_out,*)"nband_k, npw_k, npwtot_k: ",nband_k, npw_k, npwtot_k 1191 1192 ! Part taken from writewf 1193 ! Note that ind_cg_mpi_to_seq is wrong in nspinor > 1 1194 npwtot=npw_k; npwso=npw_k*nspinor 1195 npwsotot=npwso 1196 nspinortot = min(2,(1+mpi_enreg%paral_spinor)*nspinor) 1197 1198 ABI_MALLOC(ind_cg_mpi_to_seq, (npwso)) 1199 if (allocated(mpi_enreg%my_kgtab)) then 1200 ikpt_this_proc=mpi_enreg%my_kpttab(ikpt) 1201 do ispinor=1,nspinor 1202 ispinor_index=ispinor 1203 if (nproc_spinor > 1) ispinor_index = mpi_enreg%me_spinor + 1 1204 ind_cg_mpi_to_seq(1+npw_k*(ispinor-1):npw_k*ispinor)=npwtot*(ispinor_index-1) & 1205 & + mpi_enreg%my_kgtab(1:npw_k,ikpt_this_proc) 1206 end do 1207 else 1208 ABI_CHECK(nproc_fft==1, "nproc_fft !=1 and my_kgtab not allocated") 1209 ind_cg_mpi_to_seq(1:npwso) = [(ipw, ipw=1,npwso)] 1210 end if 1211 1212 ! TODO: Blocking + collective IO in comm_cell 1213 !band_blocksize = nband_k/nbdblock 1214 !ibandmin = 1 1215 !step=min(ii,MAXBAND, nband_disk) 1216 !do iblock=1,nband_disk/step+1 1217 ! ibandmax = min(ibandmin+step-1, nband_disk) 1218 ! nband_block = ibandmax-ibandmin+1 1219 !gbase = icg + npw_k*my_nspinor*nband_k 1220 1221 ! Each MPI proc write bcount bands with npwtot_k G-vectors starting from bstart. 1222 call cg2seqblocks(npwtot_k,npw_k,nband_k,cg(:,icg+1:),ind_cg_mpi_to_seq,comm_bandfft,bstart,bcount,my_cgblock) 1223 1224 ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 1225 ncerr = nf90_put_var(ncid, cg_varid, my_cgblock, start=[1,1,1,bstart,ikpt,spin], count=[2,npwtot_k,nspinor,bcount,1,1]) 1226 NCF_CHECK_MSG(ncerr, "put_var cg") 1227 1228 ABI_FREE(my_cgblock) 1229 ABI_FREE(ind_cg_mpi_to_seq) 1230 1231 icg = icg+npw_k*my_nspinor*nband_k 1232 end do !ikpt 1233 end do !spin 1234 end if 1235 1236 NCF_CHECK(nf90_close(ncid)) 1237 done = .True. 1238 1239 call cwtime_report(" scattered ncwrite", cpu, wall, gflops) 1240 #endif 1241 ! HAVE_NETCDF 1242 end if !nctk_has_mpiio 1243 end if 1244 1245 call crystal%free() 1246 if (response == 0) call ebands_free(gs_ebands) 1247 1248 ABI_FREE(occ3d) 1249 1250 DBG_EXIT("COLL") 1251 1252 end subroutine cg_ncwrite
m_iowf/kg2seqblocks [ Functions ]
[ Top ] [ m_iowf ] [ Functions ]
NAME
kg2seqblocks
FUNCTION
INPUTS
npwtot_k npw_k= kg(3,npw_k)=reduced planewave coordinates. gmpi2seq comm_fft=FFT communicator
OUTPUT
start_pwblock count_pwblock gblock(:,:)
SOURCE
1491 subroutine kg2seqblocks(npwtot_k,npw_k,kg_k,gmpi2seq,comm_fft,start_pwblock,count_pwblock,gblock) 1492 1493 !Arguments ------------------------------------ 1494 !scalars 1495 integer,intent(in) :: npwtot_k,npw_k,comm_fft 1496 integer,intent(out) :: start_pwblock,count_pwblock 1497 !arrays 1498 integer,intent(in) :: kg_k(3,npw_k),gmpi2seq(npw_k) 1499 integer,allocatable,intent(out) :: gblock(:,:) 1500 1501 !Local variables------------------------------- 1502 !scalars 1503 integer :: me_fft,nproc_fft,rank,ig,igseq,maxnpw,ierr 1504 !arrays 1505 integer,allocatable :: igstart_rank(:),gbuf(:,:) 1506 1507 ! ************************************************************************* 1508 1509 me_fft = xmpi_comm_rank(comm_fft); nproc_fft = xmpi_comm_size(comm_fft) 1510 1511 ! Handle sequential case. 1512 if (nproc_fft == 1) then 1513 start_pwblock = 1; count_pwblock = npw_k 1514 ABI_MALLOC(gblock, (3, npw_k)) 1515 gblock(:,:) = kg_k 1516 return ! DOH 1517 end if 1518 1519 ABI_MALLOC(igstart_rank, (0:nproc_fft)) 1520 igstart_rank = [(1 + (npwtot_k/nproc_fft)*rank, rank=0,nproc_fft-1), 1 + npwtot_k] 1521 1522 ! Get max dimension for workspace array 1523 ! Cannot use npwtot_k / nproc_fft because G-vectors are not equally distributed. 1524 call xmpi_max(npw_k, maxnpw, comm_fft, ierr) 1525 ABI_MALLOC(gbuf, (3, maxnpw)) 1526 1527 do rank=0,nproc_fft-1 1528 start_pwblock = igstart_rank(rank) 1529 count_pwblock = igstart_rank(rank+1) - igstart_rank(rank) 1530 1531 gbuf = 0 1532 do ig=1,npw_k 1533 igseq = gmpi2seq(ig) 1534 if (igseq >= igstart_rank(rank) .and. igseq < igstart_rank(rank+1)) then 1535 igseq = igseq - start_pwblock + 1 1536 !ABI_CHECK(igseq <= maxnpw, "boom") 1537 gbuf(:,igseq) = kg_k(:,ig) 1538 end if 1539 end do 1540 call xmpi_sum_master(gbuf,rank,comm_fft,ierr) 1541 1542 if (me_fft == rank) then 1543 ABI_MALLOC_OR_DIE(gblock, (3, count_pwblock), ierr) 1544 gblock = gbuf(:, :count_pwblock) 1545 end if 1546 end do 1547 1548 start_pwblock = igstart_rank(me_fft) 1549 count_pwblock = igstart_rank(me_fft+1) - igstart_rank(me_fft) 1550 1551 ABI_FREE(gbuf) 1552 ABI_FREE(igstart_rank) 1553 1554 end subroutine kg2seqblocks
m_iowf/ncwrite_eigen1_occ [ Functions ]
[ Top ] [ m_iowf ] [ Functions ]
NAME
ncwrite_eigen1_occ
FUNCTION
Write the first order DFPT eigenvalues and the occupations.
INPUTS
ncid=Netcdf file handler.
nband(nkpt*nsppol)=Number of bands.
mband=maximum number of bands
nkpt=Total number of k points
nsppol=1 for unpolarized, 2 for spin-polarized
eigen((2*mband**2*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point
occ3d(mband*nkpt*nsppol)=occupations for all bands at each k point
Note that occ3d differes from the occ arrays used in the rest of the code.
occ3d is an array with constant stride `mband` whereas abinit (for reasons that are not clear to me)
packs the occupations in a 1d vector with a k-dependent separation (nband_k).
occ and occ3d differ only if nband is k-dependent but you should never assume this, hence
remember to *convert* occ into occ3d before calling this routine.
SOURCE
1404 subroutine ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d) 1405 1406 !Arguments ------------------------------------ 1407 !scalars 1408 integer,intent(in) :: ncid,mband,nkpt,nsppol 1409 !arrays 1410 integer, intent(in) :: nband(nkpt*nsppol) 1411 real(dp),intent(in) :: eigen(2*mband**2*nkpt*nsppol),occ3d(mband,nkpt,nsppol) 1412 1413 !Local variables------------------------------- 1414 !scalars 1415 integer :: idx,spin,ikpt,nband_k,ib2,ib1 1416 integer :: ncerr,occ_varid,h1mat_varid,eigens_varid 1417 !arrays 1418 real(dp),allocatable :: h1mat(:,:,:,:,:), fake_eigens(:,:,:) 1419 1420 ! ************************************************************************* 1421 1422 ! Declare h1 array with abinit conventions 1423 ! Cannot use a 3D array since eigen1 are packed and one could have different number of bands 1424 ! per k-points. (this is not an official etsf-io variable)! 1425 ! elements in eigen are packed in the first positions. 1426 ! Remember that eigen are not MPI-distributed so no communication is needed 1427 idx=1 1428 ABI_CALLOC(h1mat, (2,mband,mband,nkpt,nsppol)) 1429 do spin=1,nsppol 1430 do ikpt=1,nkpt 1431 nband_k = nband(ikpt + (spin-1)*nkpt) 1432 do ib2=1,nband_k 1433 do ib1=1,nband_k 1434 h1mat(:,ib1,ib2,ikpt,spin) = eigen(idx:idx+1) 1435 idx=idx+2 1436 end do 1437 end do 1438 end do 1439 end do 1440 1441 NCF_CHECK(nctk_set_defmode(ncid)) 1442 1443 ncerr = nctk_def_arrays(ncid, nctkarr_t('h1_matrix_elements', "dp", & 1444 &"complex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")) 1445 NCF_CHECK(ncerr) 1446 1447 NCF_CHECK(nf90_inq_varid(ncid, "occupations", occ_varid)) 1448 NCF_CHECK(nf90_inq_varid(ncid, "h1_matrix_elements", h1mat_varid)) 1449 NCF_CHECK(nf90_inq_varid(ncid, "eigenvalues", eigens_varid)) 1450 1451 ! Write data 1452 NCF_CHECK(nctk_set_datamode(ncid)) 1453 NCF_CHECK_MSG(nf90_put_var(ncid, occ_varid, occ3d), "putting occ3d") 1454 NCF_CHECK_MSG(nf90_put_var(ncid, h1mat_varid, h1mat), "putting h1mat") 1455 1456 ABI_FREE(h1mat) 1457 1458 ! GS eigenvalues are set to zero. 1459 ABI_CALLOC(fake_eigens, (mband,nkpt,nsppol)) 1460 NCF_CHECK_MSG(nf90_put_var(ncid, eigens_varid, fake_eigens), "putting fake eigens") 1461 ABI_FREE(fake_eigens) 1462 1463 end subroutine ncwrite_eigen1_occ
m_iowf/outresid [ Functions ]
[ Top ] [ m_iowf ] [ Functions ]
NAME
outresid
FUNCTION
- Compute the maximal residual and eventually print it
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset kptns(3,nkpt)=k points in terms of recip primitive translations mband=maximum number of bands nband=number of bands nkpt=number of k points nsppol=1 for unpolarized, 2 for spin-polarized resid(mband*nkpt*nsppol)=squared residuals for each band and k point where resid(n,k)=|<C(n,k)|(H-e(n,k))|C(n,k)>|^2 for the ground state
OUTPUT
(only writing)
SOURCE
88 subroutine outresid(dtset,kptns,mband,nband,nkpt,nsppol,resid) 89 90 !Arguments ------------------------------------ 91 !scalars 92 integer,intent(in) :: mband,nkpt,nsppol 93 type(dataset_type),intent(in) :: dtset 94 !arrays 95 integer, intent(in) :: nband(nkpt*nsppol) 96 real(dp), intent(in) :: kptns(3,nkpt) 97 real(dp), intent(in) :: resid(mband*nkpt*nsppol) 98 99 !Local variables------------------------------- 100 integer,parameter :: nkpt_max=50 101 integer :: band_index,spin,ikpt,ibdkpt,nkpt_eff,nband_k,nband_eff 102 integer :: iband, ii 103 real(dp) :: resims,residm, residk 104 character(len=500) :: msg 105 106 !Compute mean square and maximum residual over all bands and k points and spins 107 !(disregard k point weights and occupation numbers here) 108 109 !Find largest residual over bands, k points, and spins, except for nbdbuf highest bands 110 !Already AVAILABLE in hdr ?! 111 ibdkpt=0 112 residm=zero 113 resims=zero 114 band_index=0 115 do spin=1,nsppol 116 do ikpt=1,nkpt 117 nband_k=nband(ikpt+(spin-1)*nkpt) 118 if (dtset%nbdbuf>0) then 119 nband_eff=max(1,nband_k-dtset%nbdbuf) 120 else 121 nband_eff=nband_k 122 end if 123 residm=max(residm,maxval(resid(ibdkpt+1:ibdkpt+nband_eff))) 124 resims=resims + sum(resid(ibdkpt+1:ibdkpt+nband_eff)) 125 ibdkpt=ibdkpt+nband_k 126 band_index=band_index + nband_eff 127 end do 128 end do 129 resims=resims/dble(band_index) 130 131 write(msg,'(a,2p,e12.4,a,e12.4)')' Mean square residual over all n,k,spin= ',resims,'; max=',residm 132 call wrtout([std_out, ab_out], msg) 133 134 band_index=0 135 nkpt_eff=nkpt 136 if( (dtset%prtvol==0 .or. dtset%prtvol==1) .and. nkpt_eff>nkpt_max ) nkpt_eff=nkpt_max 137 138 !Loop over spin again 139 do spin=1,nsppol 140 ! Give (squared) residuals for all bands at each k 141 do ikpt=1,nkpt 142 nband_k=nband(ikpt+(spin-1)*nkpt) 143 ! Will not print all residuals when prtvol=0 or 1 144 if(ikpt<=nkpt_eff)then 145 ! Find largest residual over all bands for given k point 146 residk=maxval(resid(1+band_index:nband_k+band_index)) 147 write(msg,'(1x,3f8.4,3x,i2,1p,e13.5,a)')kptns(1:3,ikpt),spin,residk,' kpt; spin; max resid(k); each band:' 148 if(dtset%prtvol>=2) call wrtout(ab_out, msg) 149 call wrtout(std_out, msg) 150 do ii=0,(nband_k-1)/8 151 write(msg,'(1x,1p,8e9.2)')(resid(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8)) 152 if(dtset%prtvol>=2) call wrtout(ab_out, msg) 153 call wrtout(std_out, msg) 154 end do 155 else if(ikpt==nkpt_eff+1)then 156 write(msg,'(2a)')' outwf : prtvol=0 or 1, do not print more k-points.',ch10 157 if(dtset%prtvol>=2) call wrtout(ab_out, msg) 158 call wrtout(std_out, msg) 159 end if 160 band_index=band_index+nband_k 161 end do 162 end do 163 164 end subroutine outresid
m_iowf/outwf [ Functions ]
[ Top ] [ m_iowf ] [ Functions ]
NAME
outwf
FUNCTION
Conduct output of a "wave-functions" file. - Compute the maximal residual - Then open a permanent file wff2 for final output of wf data - Create a new header for the file. - Write wave-functions (and energies)
INPUTS
cg(2,mcg)=wavefunction array (storage if nkpt>1) dtset <type(dataset_type)>=all input variables for this dataset psps <type(pseudopotential_type)>=variables related to pseudopotentials eigen( (2*mband)**response *mband*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point filnam= character string giving the root to form the name of the output WFK or WFQ file if response==0, otherwise it is the filename. hdr <type(hdr_type)>=the header of wf, den and pot files kg(3,mpw*mkmem)=reduced planewave coordinates. kptns(3,nkpt)=k points in terms of recip primitive translations mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mkmem=Number of k-points treated by this node. mpi_enreg=information about MPI parallelization mpw=maximum number of plane waves natom=number of atoms in unit cell nband=number of bands nkpt=number of k points npwarr(nkpt)=number of planewaves in basis and on boundary for each k nsppol=1 for unpolarized, 2 for spin-polarized nstep=desired number of electron iteration steps occ(mband*nkpt*nsppol)=occupations for all bands at each k point response: if == 0, GS wavefunctions , if == 1, RF wavefunctions unwff2=unit for output of wavefunction wfs <type(wvl_projector_type)>=wavefunctions information for wavelets.
OUTPUT
(only writing)
NOTES
* The name of the file wff2 might be the same as that of the file wff1.
SOURCE
212 subroutine outwf(cg,dtset,psps,eigen,filnam,hdr,kg,kptns,mband,mcg,mkmem,& 213 mpi_enreg,mpw,natom,nband,nkpt,npwarr,& 214 nsppol,occ,response,unwff2,& 215 wfs,wvl,force_write) ! optional 216 217 !Arguments ------------------------------------ 218 !scalars 219 integer,intent(in) :: mband,mcg,mkmem,mpw,natom,nkpt,nsppol,response,unwff2 220 !integer,intent(in) :: nstep 221 character(len=*),intent(in) :: filnam 222 type(MPI_type),intent(in) :: mpi_enreg 223 type(dataset_type),intent(in) :: dtset 224 type(pseudopotential_type),intent(in) :: psps 225 type(hdr_type), intent(inout) :: hdr 226 type(wvl_wf_type),intent(in) :: wfs 227 type(wvl_internal_type), intent(in) :: wvl 228 logical, intent(in), optional :: force_write 229 !arrays 230 integer, intent(in) :: kg(3,mpw*mkmem),nband(nkpt*nsppol),npwarr(nkpt) 231 real(dp), intent(inout) :: cg(2,mcg) 232 real(dp), intent(in) :: eigen((2*mband)**response*mband*nkpt*nsppol),kptns(3,nkpt) 233 real(dp), intent(in) :: occ(mband*nkpt*nsppol) 234 235 !Local variables------------------------------- 236 integer :: iomode,action,band_index,fform,formeig,iband,icg !,iat,iproj 237 integer :: ierr,ikg,ikpt,spin,master,mcg_disk,me,me0,mtag,my_nspinor 238 integer :: nband_k,nmaster,npw_k,option,rdwr,sender,source !npwtot_k, 239 integer :: spaceComm,spaceComm_io,spacecomsender,spaceWorld,sread,sskip,tim_rwwf,xfdim2 240 #ifdef HAVE_MPI 241 integer :: ipwnbd 242 #endif 243 real(dp) :: cpu,wall,gflops 244 logical :: ihave_data,iwrite,iam_master,done,prtwf 245 character(len=500) :: msg 246 type(wffile_type) :: wff2 247 !character(len=fnlen) :: path 248 !arrays 249 integer,allocatable :: kg_disk(:,:) 250 real(dp) :: tsec(2) 251 real(dp),allocatable :: cg_disk(:,:),eig_k(:),occ_k(:) 252 253 ! ************************************************************************* 254 !For readability of the source file, define a "me" variable also in the sequential case 255 256 DBG_ENTER("COLL") 257 258 ABI_UNUSED(kptns(1,1)) 259 260 xfdim2 = natom+4 261 !Init mpi_comm 262 spaceWorld= mpi_enreg%comm_cell 263 spaceComm=spaceWorld 264 spaceComm_io=xmpi_comm_self 265 266 if (mpi_enreg%paral_kgb==1 ) spaceComm_io= mpi_enreg%comm_bandspinorfft 267 if (mpi_enreg%paral_kgb==1 ) spaceComm= mpi_enreg%comm_cell 268 269 !Paral_kgb=1 and Fortran-I/O is not supported (only for testing purpose) 270 if (mpi_enreg%paral_kgb==1.and.dtset%iomode==IO_MODE_FORTRAN) then 271 spaceWorld=mpi_enreg%comm_kpt 272 write(msg,'(7a)') & 273 'WF file is written using standard Fortran I/O',ch10,& 274 'and Kpt-band-FFT parallelization is active !',ch10,& 275 'This is only allowed for testing purposes.',ch10,& 276 'The produced WF file will be incomplete and not useable.' 277 ABI_WARNING(msg) 278 end if 279 280 !If parallel HF calculation 281 if (mpi_enreg%paral_hf==1 ) spaceComm_io= mpi_enreg%comm_hf 282 if (mpi_enreg%paral_hf==1 ) spaceComm= mpi_enreg%comm_cell 283 284 !Paral_hf=1 and Fortran-I/O is not supported (copy from paral_kgb... not tested) 285 if (mpi_enreg%paral_hf==1.and.dtset%iomode==IO_MODE_FORTRAN) then 286 spaceWorld=mpi_enreg%comm_kpt 287 write(msg,'(7a)') & 288 'WF file is written using standard Fortran I/O',ch10,& 289 'and HF parallelization is active !',ch10,& 290 'This is only allowed for testing purposes.',ch10,& 291 'The produced WF file will be incomplete and not useable.' 292 ABI_WARNING(msg) 293 end if 294 295 ! check consistency between dimensions and input hdr. 296 !ABI_CHECK(mband == maxval(hdr%nband), "hdr:mband") 297 !ABI_CHECK(nkpt == hdr%nkpt, "hdr:nkpt") 298 !ABI_CHECK(nsppol == hdr%nsppol, "hdr:nsppol") 299 !ABI_CHECK(all(hdr%npwarr == npwarr), "hdr:npwarr") 300 !ABI_CHECK(all(hdr%nband == nband), "hdr:nband") 301 !ABI_CHECK(maxval(hdr%npwarr) == mpw, "hdr:nband") 302 303 !Init me 304 me=mpi_enreg%me_kpt 305 me0=me 306 !Define master 307 master=0 308 309 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 310 tim_rwwf =0 311 source = master 312 sread = master 313 iam_master=(master==me) 314 iwrite=iam_master 315 sender=-1 316 317 !Will write the wavefunction file only when nstep>0 318 !MT 07 2015: writing reactivated when nstep=0 319 !if (nstep>0 .and. dtset%prtwf/=0) then 320 !FB 03/2022: Added an option to force writing (used in RT-TDDFT) 321 prtwf = dtset%prtwf/=0 322 if (present(force_write)) then 323 if (force_write) prtwf = .true. 324 end if 325 if (prtwf) then 326 327 ! Only the master write the file, except if MPI I/O, but the 328 ! full wff dataset should be provided to WffOpen in this case 329 iomode=IO_MODE_FORTRAN_MASTER 330 if (dtset%iomode==IO_MODE_MPI) iomode = IO_MODE_MPI 331 if (dtset%iomode==IO_MODE_ETSF) iomode = IO_MODE_ETSF 332 333 #ifdef HAVE_NETCDF 334 if (dtset%iomode == IO_MODE_ETSF .and. dtset%usewvl == 0) then 335 call cg_ncwrite(filnam,hdr,dtset,response,mpw,mband,nband,nkpt,nsppol,& 336 dtset%nspinor,mcg,mkmem,eigen,occ,cg,npwarr,kg,mpi_enreg,done) 337 #ifdef HAVE_NETCDF_DEFAULT 338 ABI_CHECK(done, "cg_ncwrite must handle the output of the WFK file.") 339 #endif 340 341 ! Write KB form factors. Only master works. G-vectors are read from file to avoid 342 ! having to deal with paral_kgb distribution. 343 if (me == master .and. dtset%prtkbff == 1 .and. dtset%iomode == IO_MODE_ETSF .and. dtset%usepaw == 0) then 344 ABI_CHECK(done, "cg_ncwrite was not able to generate WFK.nc in parallel. Perhaps hdf5 is not working") 345 call prtkbff(filnam, hdr, psps, dtset%prtvol) 346 end if 347 348 if (done) return 349 ! If cg_ncwrite cannot handle the IO because HDF5 + MPI-IO support is missing, we fallback to Fortran + MPI-IO. 350 msg = "Could not produce a netcdf file in parallel (MPI-IO support is missing). Will fallback to MPI-IO with Fortran" 351 ABI_WARNING(msg) 352 iomode=IO_MODE_MPI 353 end if 354 #endif 355 356 call cwtime(cpu, wall, gflops, "start") 357 call wrtout(std_out, sjoin(ch10,' outwf: writing wavefunctions to:', trim(filnam), "with iomode:", iomode2str(iomode))) 358 359 ! Create an ETSF file for the wavefunctions 360 if (iomode == IO_MODE_ETSF) then 361 ABI_CHECK(xmpi_comm_size(spaceComm) == 1, "Legacy etsf-io code does not support nprocs > 1") 362 ABI_ERROR("ETSF_IO is not activated") 363 ABI_UNUSED(psps%ntypat) 364 end if 365 366 call WffOpen(iomode,spaceComm,filnam,ierr,wff2,master,me0,unwff2,spaceComm_io) 367 ! Conduct wavefunction output to wff2 368 369 ABI_MALLOC(kg_disk,(3,mpw)) 370 371 mcg_disk=mpw*my_nspinor*mband 372 formeig=0; if (response==1) formeig=1 373 374 ABI_MALLOC(eig_k,( (2*mband)**formeig * mband)) 375 ABI_MALLOC(occ_k,(mband)) 376 377 #ifdef HAVE_MPI 378 call xmpi_barrier(spaceComm) 379 ! Compute mband and mpw 380 ABI_MALLOC_OR_DIE(cg_disk,(2,mcg_disk), ierr) 381 #endif 382 383 band_index=0 384 icg=0 385 if(mpi_enreg%paralbd==0) tim_rwwf=6 386 if(mpi_enreg%paralbd==1) tim_rwwf=12 387 388 ! Write header info for new wf file 389 rdwr=2 390 if (dtset%usewvl==0) then 391 fform=2 392 else 393 fform = 200 ! Use 200 as radical for naming file format used by wavelets. 394 end if 395 396 if (wff2%iomode < 2) then 397 call hdr_io(fform,hdr,rdwr,wff2) 398 call WffKg(wff2,1) 399 else if (wff2%iomode==IO_MODE_ETSF .and. iam_master) then 400 #ifdef HAVE_NETCDF 401 NCF_CHECK(hdr%ncwrite(wff2%unwff, fform, nc_define=.True.)) 402 #endif 403 end if 404 405 do spin=1,nsppol 406 ikg=0 407 408 do ikpt=1,nkpt 409 nband_k=nband(ikpt+(spin-1)*nkpt) 410 npw_k=npwarr(ikpt) 411 412 #ifdef HAVE_MPI 413 if (dtset%usewvl == 0) then 414 mtag=ikpt+(spin-1)*nkpt 415 call xmpi_barrier(spaceWorld) 416 417 ! Must transfer the wavefunctions to the master processor 418 ! Separate sections for paralbd=1 or other values ; might be merged 419 if(mpi_enreg%paralbd==0)then 420 nmaster=0 421 source=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,spin)) 422 ihave_data=.false. 423 if(source==me)ihave_data=.true. 424 action=0 425 ! I am the master node, and I have the data in cg or cg_disk 426 if((iam_master).and.(ihave_data))action=1 427 ! I am not the master, and I have the data => send to master 428 if((.not.iam_master).and.(ihave_data))action=2 429 ! I am the master, and I receive the data 430 if((iam_master).and.(.not.ihave_data))action=3 431 432 ! I have the data in cg or cg_disk ( MPI_IO case) 433 if (iomode==IO_MODE_MPI) then 434 action = 0 435 sender=-1 436 iwrite=.false. 437 if (ihave_data)then 438 action=1 439 iwrite=.true. 440 sender=me 441 end if 442 end if 443 444 ! I am the master node, and I have the data in cg or cg_disk 445 ! I have the data in cg or cg_disk ( MPI_IO case) 446 if(action==1)then 447 ! Copy from kg to kg_disk 448 kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 449 ! Copy from cg to cg_disk 450 do ipwnbd=1,nband_k*npw_k*my_nspinor 451 cg_disk(1,ipwnbd)=cg(1,ipwnbd+icg) 452 cg_disk(2,ipwnbd)=cg(2,ipwnbd+icg) 453 end do 454 end if 455 456 ! I am not the master, and I have the data => send to master 457 ! I am the master, and I receive the data 458 if ( action==2.or.action==3) then 459 !write(std_out,*)npw_k,nband_k 460 call timab(48,1,tsec) 461 if(action==2)then 462 call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,source,kg_disk,nmaster,spaceWorld,2*mtag+1,ierr) 463 call xmpi_exch(cg(:,icg+1:icg+nband_k*npw_k*my_nspinor),2*nband_k*npw_k*my_nspinor, & 464 & source,cg_disk,nmaster,spaceWorld,2*mtag+2,ierr) 465 else 466 call xmpi_exch(kg_disk,3*npw_k,source,kg_disk,nmaster,spaceWorld,2*mtag+1,ierr) 467 call xmpi_exch(cg_disk,2*nband_k*npw_k*my_nspinor,source,cg_disk,nmaster,spaceWorld,2*mtag+2,ierr) 468 end if 469 call timab(48,2,tsec) 470 end if 471 472 473 else if(mpi_enreg%paralbd==1)then 474 nmaster=0 475 #ifdef HAVE_MPI_IO 476 sender=IO_MODE_FORTRAN_MASTER 477 if( iomode==IO_MODE_MPI) then 478 nmaster=mpi_enreg%proc_distrb(ikpt,1,spin) 479 sender=nmaster 480 end if 481 #endif 482 483 ! Note the loop over bands 484 do iband=1,nband_k 485 486 ! The message passing related to kg is counted as one band 487 action=0 488 489 ! I am the master node, and I have the data in cg or cg_disk 490 if( mpi_enreg%proc_distrb(ikpt,iband,spin)==nmaster .and. me==nmaster) then 491 action=1 492 ! I am not the master, and I have the data => send to master 493 elseif( mpi_enreg%proc_distrb(ikpt,iband,spin)==me .and. me/=nmaster ) then 494 action = 2 495 ! I am the master, and I receive the data 496 elseif( mpi_enreg%proc_distrb(ikpt,iband,spin)/=me .and. me==nmaster ) then 497 action=3 498 end if 499 500 if(action==1) then 501 ! I am the master node, and I have the data in cg or cg_disk 502 ! Copy from kg to kg_disk 503 if(iband==1)kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 504 ! Copy from cg to cg_disk 505 do ipwnbd=1,npw_k*my_nspinor 506 cg_disk(1,(iband-1)*npw_k*my_nspinor+ipwnbd) = cg(1,(iband-1)*npw_k*my_nspinor+ipwnbd+icg) 507 cg_disk(2,(iband-1)*npw_k*my_nspinor+ipwnbd) = cg(2,(iband-1)*npw_k*my_nspinor+ipwnbd+icg) 508 end do 509 end if ! action=1 510 511 if ( action==2.or.action==3) then 512 ! action=2 : I am not the master, and I have the data => send to master 513 ! action=3 : I am the master, and I receive the data 514 call timab(48,1,tsec) 515 if ( iband == 1 ) then 516 if (action==2) then 517 call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,mpi_enreg%proc_distrb(ikpt,iband,spin), & 518 & kg_disk,nmaster,spaceWorld,iband*(mtag-1)+1,ierr) 519 else 520 call xmpi_exch(kg_disk,3*npw_k,mpi_enreg%proc_distrb(ikpt,iband,spin), & 521 & kg_disk,nmaster,spaceWorld,iband*(mtag-1)+1,ierr) 522 end if 523 end if ! iband =1 524 ipwnbd=(iband-1)*npw_k*my_nspinor 525 if (action==2) then 526 call xmpi_exch( cg(:,ipwnbd+icg+1:ipwnbd+icg+npw_k*my_nspinor),2*npw_k*my_nspinor & 527 & ,mpi_enreg%proc_distrb(ikpt,iband,spin) & 528 & ,cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),nmaster,spaceWorld,iband*(mtag-1)+2,ierr) 529 else 530 call xmpi_exch( cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),2*npw_k*my_nspinor & 531 & ,mpi_enreg%proc_distrb(ikpt,iband,spin) & 532 & ,cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),nmaster,spaceWorld,iband*(mtag-1)+2,ierr) 533 end if 534 535 call timab(48,2,tsec) 536 end if ! action=2 or action=3 537 538 if(iomode==IO_MODE_MPI) then 539 ! I have the data in cg or cg_disk 540 iwrite=.false. 541 if (nmaster == me) iwrite=.true. 542 end if 543 544 end do ! End of loop over bands 545 end if ! End of paralbd=1 546 end if 547 #endif 548 549 ! Only the master will write to disk the final output wf file. 550 ! in MPI_IO case only iwrite will write to disk the final output wf file. 551 if(iwrite) then 552 ! write(std_out,*) 'outwf : I am master and will write wf file' 553 if(formeig==0)then 554 eig_k(1:nband_k)=eigen(1+band_index:nband_k+band_index) 555 occ_k(1:nband_k)=occ(1+band_index:nband_k+band_index) 556 else 557 eig_k(1:2*nband_k*nband_k)=eigen(1+band_index:2*nband_k*nband_k+band_index) 558 end if 559 option=2 560 if(dtset%prtwf==3)option=5 561 ! if (dtset%prtwf == 2 .and. mkmem/=0) option=4 562 563 if (dtset%usewvl == 0) then 564 #ifdef HAVE_MPI 565 call rwwf(cg_disk,eig_k,formeig,0,0,ikpt,spin,kg_disk,mband,mcg_disk,mpi_enreg, & 566 & nband_k, nband_k,npw_k,my_nspinor,occ_k,option,1,tim_rwwf,wff2) 567 568 #else 569 kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 570 call rwwf(cg,eig_k,formeig,0,icg,ikpt,spin,kg_disk,mband,mcg,mpi_enreg,nband_k, & 571 & nband_k, npw_k,my_nspinor,occ_k,option,1,tim_rwwf,wff2) 572 #endif 573 else 574 call wvl_write(dtset,eigen,mpi_enreg,option,hdr%rprimd,wff2,wfs,wvl,hdr%xred) 575 end if 576 end if 577 578 ! The wavefunctions for the present k point and spin are written 579 if(response==0)band_index=band_index+nband_k 580 if(response==1)band_index=band_index+2*nband_k*nband_k 581 582 sskip=1 583 #ifdef HAVE_MPI 584 if (dtset%usewvl == 0) then 585 sskip=0 586 if(.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me)))sskip=1 587 end if 588 #endif 589 if(sskip==1)then 590 icg=icg+npw_k*my_nspinor*nband_k 591 ikg=ikg+npw_k 592 end if 593 594 595 #ifdef HAVE_MPI_IO 596 spacecomsender=spaceComm 597 if (mpi_enreg%paral_kgb==1) spacecomsender =mpi_enreg%comm_kpt 598 if (mpi_enreg%paral_hf==1) spacecomsender =mpi_enreg%comm_kpt 599 call WffOffset(wff2,sender,spacecomsender,ierr) 600 #endif 601 602 end do ! ikpt 603 end do ! spin 604 605 ABI_FREE(kg_disk) 606 #ifdef HAVE_MPI 607 ABI_FREE(cg_disk) 608 #endif 609 610 ABI_FREE(eig_k) 611 ABI_FREE(occ_k) 612 613 ! Close the wavefunction file (and do NOT delete it !) 614 !if (wff2%iomode /= IO_MODE_NETCDF) then 615 call WffClose(wff2,ierr) 616 !end if 617 618 call cwtime_report(" WFK output", cpu, wall, gflops) 619 end if ! End condition of nstep>0 620 621 ! Block here because we might need to read the WFK file in the caller. 622 call xmpi_barrier(mpi_enreg%comm_cell) 623 624 DBG_EXIT("COLL") 625 626 end subroutine outwf
m_iowf/prtkbff [ Functions ]
[ Top ] [ m_iowf ] [ Functions ]
NAME
prtkbff
FUNCTION
Write KB form factors to WFK in netcdf format. Only master works. G-vectors are read from file to avoid having to deal with paral_kgb distribution.
INPUTS
OUTPUT
SOURCE
1270 subroutine prtkbff(filnam, hdr, psps, prtvol) 1271 1272 !Arguments ------------------------------------ 1273 character(len=*),intent(in) :: filnam 1274 type(hdr_type),intent(in) :: hdr 1275 type(pseudopotential_type),intent(in) :: psps 1276 integer,intent(in) :: prtvol 1277 1278 !Local variables------------------------------- 1279 !scalars 1280 character(len=fnlen) :: path 1281 character(len=500) :: msg 1282 integer :: ncid, ncerr, kg_varid, mpw_disk, npwk_disk, ikpt, iat, iproj 1283 integer,allocatable :: kg_disk(:,:) 1284 real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:),vkbsign(:,:) 1285 type(crystal_t) :: crystal 1286 1287 ! ************************************************************************* 1288 1289 path = nctk_ncify(filnam) 1290 call wrtout(std_out, sjoin("Writing KB form factors to:", path)) 1291 NCF_CHECK(nctk_open_modify(ncid, path, xmpi_comm_self)) 1292 NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid)) 1293 mpw_disk = maxval(hdr%npwarr) 1294 1295 ! Dimensions needed by client code to allocate memory when reading. 1296 ncerr = nctk_def_dims(ncid, [ & 1297 nctkdim_t("mproj", psps%mproj), & 1298 nctkdim_t("mpsang", psps%mpsang), & 1299 nctkdim_t("mpssoang", psps%mpssoang), & 1300 nctkdim_t("lnmax", psps%lnmax), & 1301 nctkdim_t("lmnmax", psps%lmnmax) & 1302 ]) 1303 NCF_CHECK(ncerr) 1304 1305 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "mpspso"]) 1306 NCF_CHECK(ncerr) 1307 1308 ! Write indlmn table (needed to access vkb arrays) 1309 ncerr = nctk_def_arrays(ncid, [ & 1310 nctkarr_t("indlmn", "int", "six, lmnmax, number_of_atom_species"), & 1311 nctkarr_t("vkbsign", "dp", "lnmax, number_of_atom_species"), & 1312 nctkarr_t("vkb", "dp", "max_number_of_coefficients, lnmax, number_of_atom_species, number_of_kpoints"), & 1313 nctkarr_t("vkbd", "dp", "max_number_of_coefficients, lnmax, number_of_atom_species, number_of_kpoints") & 1314 ], defmode=.True.) 1315 NCF_CHECK(ncerr) 1316 1317 ! Switch to write mode. 1318 NCF_CHECK(nctk_set_datamode(ncid)) 1319 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "indlmn"), psps%indlmn)) 1320 1321 ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: & 1322 "mpspso"], & 1323 [psps%mpspso]) 1324 NCF_CHECK(ncerr) 1325 1326 ! Calculate KB form factors and derivatives. 1327 ! The arrays are allocated with lnmax to support pseudos with more than projector. 1328 ! Note that lnmax takes into account lloc hence arrays are in packed form and should be 1329 ! accessed with the indices provided by psps%indlmn. 1330 ABI_MALLOC(vkbsign, (psps%lnmax, psps%ntypat)) 1331 ABI_MALLOC(vkb, (mpw_disk, psps%lnmax, psps%ntypat)) 1332 ABI_MALLOC(vkbd, (mpw_disk, psps%lnmax, psps%ntypat)) 1333 ABI_MALLOC(kg_disk, (3, mpw_disk)) 1334 1335 crystal = hdr%get_crystal() 1336 1337 ! For each k-point: read full G-vector list from file, compute KB data and write to file. 1338 do ikpt=1,hdr%nkpt 1339 npwk_disk = hdr%npwarr(ikpt) 1340 NCF_CHECK(nf90_get_var(ncid, kg_varid, kg_disk, start=[1, 1, ikpt], count=[3, npwk_disk, 1])) 1341 vkb = zero; vkbd = zero 1342 call calc_vkb(crystal, psps, hdr%kptns(:, ikpt), npwk_disk, mpw_disk, kg_disk, vkbsign, vkb, vkbd) 1343 1344 if (ikpt == 1) then 1345 ! This for the automatic tests. 1346 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkbsign"), vkbsign)) 1347 if(prtvol >= 2) then 1348 write(msg,'(a)') 'prtkbff: writing first and last G-components of the KB form factors' 1349 call wrtout(ab_out, msg) 1350 do iat=1,psps%ntypat 1351 write(msg,'(a10,i5)') 'atype ',iat 1352 call wrtout(ab_out, msg) 1353 do iproj=1,psps%lnmax 1354 write(msg,'(a10,i5,a,a10,e12.4,a,2(a10,2e12.4,a))') & 1355 'projector ', iproj,ch10, & 1356 'vkbsign ', vkbsign(iproj,iat), ch10, & 1357 'vkb ', vkb(1,iproj,iat), vkb(npwk_disk,iproj,iat), ch10, & 1358 'vkbd ', vkbd(1,iproj,iat), vkbd(npwk_disk,iproj,iat), '' 1359 call wrtout(ab_out, msg) 1360 end do 1361 end do 1362 end if 1363 end if 1364 1365 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkb"), vkb, start=[1, 1, 1, ikpt])) 1366 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkbd"), vkbd, start=[1, 1, 1, ikpt])) 1367 end do 1368 NCF_CHECK(nf90_close(ncid)) 1369 1370 ABI_FREE(kg_disk) 1371 ABI_FREE(vkbsign) 1372 ABI_FREE(vkb) 1373 ABI_FREE(vkbd) 1374 call crystal%free() 1375 1376 end subroutine prtkbff