TABLE OF CONTENTS
- ABINIT/dsksta
- ABINIT/m_io_kss
- ABINIT/memkss
- m_io_kss/k2gamma_centered
- m_io_kss/kss_calc_vkb
- m_io_kss/make_gvec_kss
- m_io_kss/outkss
- m_io_kss/write_kss_header
- m_io_kss/write_kss_wfgk
- m_io_kss/write_vkb
ABINIT/dsksta [ Functions ]
NAME
dsksta
FUNCTION
This routine evaluates the amount of disk space required by the _KSS file.
INPUTS
dimlmn(natom*usepaw)=Number of nlm partial waves for each atom. ishm=Number of G-shells to be saved in _KSS file. mpsang=Max angular momentum +1 for pseudos. natom=Number of atoms in the unit cell. nbandkss=Number of desired bands to be saved in _KSS file nkpt=Number of k points. npwkss=Number of desired G-vectors to be saved in _KSS file. nspinor=Number of spinorial components. nsppol=Number of independent spin polarizations. ntypat=Number of type of atoms. nsym2=Number of symmetries in space group, without INV usepaw=1 if PAW.
OUTPUT
Writes on standard output
SOURCE
2118 subroutine dsksta(ishm,usepaw,nbandkss,mpsang,natom,ntypat,npwkss,nkpt,nspinor,nsppol,nsym2,dimlmn) 2119 2120 !Arguments ------------------------------------ 2121 !scalars 2122 integer,intent(in) :: usepaw,ishm,nbandkss,mpsang,natom,ntypat,nkpt 2123 integer,intent(in) :: npwkss,nspinor,nsppol,nsym2 2124 !arrays 2125 integer,intent(in) :: dimlmn(natom*usepaw) 2126 2127 !Local variables------------------------------- 2128 !scalars 2129 integer :: bsize_tot,bsize_hdr,bsize_kb,bsize_wf,bsize_cprj 2130 character(len=500) :: msg 2131 2132 ! ********************************************************************* 2133 2134 !The Abinit header is not considered. 2135 bsize_hdr= 80*2 + & !title 2136 &5*4 + & !nsym2,nbandksseff,npwkss,ishm,mpsang 2137 &nsym2*9*4 + & !symrel2 2138 &nsym2*3*8 + & !tnons 2139 &npwkss*3*4 + & !gbig 2140 &ishm*4 !shlim 2141 2142 !NOTE: vkb does not depend on nsppol, however the elements are written for each spin. 2143 bsize_kb=0 2144 if (usepaw==0) then 2145 bsize_kb= nsppol* & 2146 & ( mpsang*ntypat *8 + & !vkbsign 2147 & 2*(nkpt*mpsang*ntypat*npwkss*8) & !vkbd,vkbd 2148 & ) 2149 end if 2150 2151 bsize_wf= nsppol* & 2152 & ( nkpt*nbandkss *8 + & !energies 2153 &nkpt*nbandkss*nspinor*npwkss*2*8 & !wfg 2154 &) 2155 2156 !For PAW add space required by projectors. 2157 bsize_cprj=0 2158 if (usepaw==1) then 2159 bsize_cprj=SUM(dimlmn(:))*(nsppol*nkpt*nspinor*nbandkss*2*8) 2160 end if 2161 2162 bsize_tot = bsize_hdr + bsize_kb + bsize_wf + bsize_cprj 2163 write(msg,'(2a,f8.2,4a,4(a,f8.2,2a))')ch10,& 2164 & ' Total amount of disk space required by _KSS file = ',bsize_tot*b2Mb,' Mb.',ch10,& 2165 & ' Subdivided into : ',ch10,& 2166 & ' Header = ',bsize_hdr *b2Mb,' Mb.',ch10,& 2167 & ' KB elements = ',bsize_kb *b2Mb,' Mb.',ch10,& 2168 & ' Wavefunctions (PW) = ',bsize_wf *b2Mb,' Mb.',ch10,& 2169 & ' PAW projectors = ',bsize_cprj*b2Mb,' Mb.',ch10 2170 call wrtout(std_out,msg,'COLL') 2171 2172 end subroutine dsksta
ABINIT/m_io_kss [ Modules ]
NAME
m_io_kss
FUNCTION
This module contains procedured dealing with the IO of the KSS file.
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (MG, MT, VO, AR, LR, RWG, MM, XG, RShaltaf) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_io_kss 24 25 use defs_basis 26 use m_abicore 27 use m_xmpi 28 use m_errors 29 use m_nctk 30 #ifdef HAVE_NETCDF 31 use netcdf 32 #endif 33 use m_hdr 34 use m_wfk 35 use m_cgtools 36 use m_hamiltonian 37 use m_electronpositron 38 use m_pawtab 39 use m_paw_ij 40 use m_pawcprj 41 use m_pawfgr 42 use m_dtfil 43 use m_dtset 44 45 use defs_datatypes, only : pseudopotential_type 46 use defs_abitypes, only : MPI_type 47 use m_time, only : timab 48 use m_io_tools, only : open_file 49 use m_fstrings, only : sjoin, itoa, strcat 50 use m_hide_lapack, only : xheevx_cplex, xhegvx_cplex 51 use m_geometry, only : metric, remove_inversion 52 use m_mpinfo, only : destroy_mpi_enreg, proc_distrb_cycle 53 use m_fftcore, only : get_kg, sphere 54 use m_fft, only : fftpac 55 use m_crystal , only : crystal_t 56 use m_gsphere, only : table_gbig2kg, merge_and_sort_kg 57 use m_kg, only : mkkin, mkkpg 58 use m_ksdiago, only : ksdiago, init_ddiago_ctl, ddiago_ctl_type 59 use m_mkffnl, only : mkffnl 60 use m_getghc, only : getghc 61 use m_fourier_interpol, only : transgrid 62 63 implicit none 64 65 private 66 67 public :: write_kss_header ! Writes the header of the KSS file. 68 !private :: write_vkb ! Writes the KB form factors and derivates on file for a single k-point. 69 public :: write_kss_wfgk ! Write the Gamma-centered wavefunctions and energies on the KSS file for a single k-point. 70 public :: k2gamma_centered ! Convert a set of wavefunctions from the k-centered to the gamma-centered basis set. 71 public :: make_gvec_kss ! Build the list of G-vectors for the KSS file. 72 public :: outkss ! Generate KSS file 73 74 CONTAINS !===========================================================
ABINIT/memkss [ Functions ]
NAME
memkss
FUNCTION
This routine evaluates the additional amount of memory required by routine 'outkss'.
INPUTS
mband=maximum number of bands mgfft=maximum size of 1D FFTs mproj=maximum dimension for number of projection operators for each angular momentum for nonlocal pseudopotential mpsang= 1+maximum angular momentum for nonlocal pseudopotentials mpw=maximum dimensioned size of npw. natom=number of atoms in cell. ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkpt=number of k points. nspinor=number of spinorial components of the wavefunctions nsym=number of symmetries in space group ntypat=number of types of atoms in unit cell.
NOTES
This routine is not available for paw calculations
SOURCE
1995 subroutine memkss(mband,mgfft,mproj,mpsang,mpw,natom,ngfft,nkpt,nspinor,nsym,ntypat) 1996 1997 !Arguments ------------------------------------ 1998 !scalars 1999 integer,intent(in) :: mband,mgfft,mproj,mpsang,mpw,natom,nkpt,nspinor 2000 integer,intent(in) :: nsym,ntypat 2001 !arrays 2002 integer,intent(in) :: ngfft(18) 2003 2004 !Local variables------------------------------- 2005 !scalars 2006 integer(i8b) :: isize,memsize 2007 character(len=500) :: msg 2008 2009 ! ********************************************************************* 2010 ! 2011 isize=580+fnlen+4*(81+nkpt+9*nsym)+8*15 !non allocatable var. 2012 if(xmpi_paral==1)then 2013 isize=isize+4*4 !kpt_distrb 2014 end if 2015 memsize=isize 2016 isize=isize+4*nkpt+12*mpw+20*nkpt*mpw !nbasek,gbasek,cnormk,gcurr 2017 memsize=max(memsize,isize) 2018 if(xmpi_paral==1)then 2019 isize=isize+12*mpw*nkpt !ibuf1,ibuf2,rbuf1 2020 memsize=max(memsize,isize) 2021 isize=isize-12*mpw*nkpt !ibuf1,ibuf2,rbuf1 2022 end if 2023 isize=isize+40*mpw !gbase,cnorm 2024 memsize=max(memsize,isize) 2025 isize=isize-4*nkpt-20*mpw*nkpt !nbasek,gbasek,cnormk 2026 isize=isize+4*mpw !insort 2027 memsize=max(memsize,isize) 2028 isize=isize-16*mpw !cnorm 2029 isize=isize+28*mpw+24*nsym !gbig,nshell,gshell 2030 memsize=max(memsize,isize) 2031 isize=isize+4*mpw !shlim 2032 memsize=max(memsize,isize) 2033 isize=isize-44*mpw-24*nsym !gcurr,gbase,gshell,insort,nshell 2034 isize=isize-4*mpw !shlim 2035 isize=isize+8*mpw*nspinor& 2036 & +16*mpw*nspinor*(mpw*nspinor+1) !eigval,eigvec 2037 memsize=max(memsize,isize) 2038 isize=isize+8*mpw+8*ngfft(4)& 2039 & *ngfft(5)*ngfft(6) !ts,vlocal 2040 memsize=max(memsize,isize) 2041 isize=isize+8*mgfft+4+28*mpw !gbound,indpw_k,kg_k 2042 memsize=max(memsize,isize) 2043 isize=isize+8*natom& 2044 & +24*mpw*ntypat*mpsang*mproj !phkxred,ffnl,kinpw 2045 memsize=max(memsize,isize) 2046 isize=isize+16*mpw*natom !ph3d 2047 memsize=max(memsize,isize) 2048 isize=isize+48*mpw*nspinor& 2049 & +8*mpw*nspinor*(mpw*nspinor+1) !pwave,subghg,gvnlg 2050 if (nspinor==2)& 2051 & isize=isize+40*mpw*nspinor !pwave_so,subghg_so 2052 memsize=max(memsize,isize) 2053 isize=isize+8*mpw*nspinor*(mpw*nspinor+1) !ghg 2054 memsize=max(memsize,isize) 2055 isize=isize+8*ngfft(4)*ngfft(5)*ngfft(6) !work 2056 memsize=max(memsize,isize) 2057 isize=isize-8*ngfft(4)*ngfft(5)*ngfft(6) !work 2058 isize=isize-8*mgfft+4+28*mpw& !gbound,indpw_k,kg_k 2059 &-8*natom-24*mpw*ntypat*mpsang*mproj& !phkxred,ffnl,kinpw 2060 &-16*mpw*natom !ph3d 2061 isize=isize-48*mpw*nspinor& 2062 & -8*mpw*nspinor*(mpw*nspinor+1) !pwave,subghg,gvnlg 2063 if (nspinor==2)& 2064 & isize=isize-40*mpw*nspinor !pwave_so,subghg_so 2065 2066 isize=isize+56*mpw*nspinor !cwork,rwork 2067 memsize=max(memsize,isize) 2068 isize=isize-56*mpw*nspinor !cwork,rwork 2069 isize=isize+112*mpw*nspinor !cwork,rwork,iwork,ifail 2070 memsize=max(memsize,isize) 2071 isize=isize-112*mpw*nspinor !cwork,rwork,iwork,ifail 2072 isize=isize-8*mpw*nspinor*(mpw*nspinor+1) !ghg 2073 isize=isize+8*mband !occ_k 2074 memsize=max(memsize,isize) 2075 isize=isize-8*mband !occ_k 2076 isize=isize-8*mpw*nspinor& 2077 & -16*mpw*nspinor*(mpw*nspinor+1) !eigval,eigvec 2078 isize=isize-32*mpw-8*ngfft(4)& 2079 & *ngfft(5)*ngfft(6) !gbig,ts,vlocal 2080 if(xmpi_paral==1)then 2081 isize=isize-4*4 !kpt_distrb 2082 end if 2083 isize=isize-580-fnlen-4*(81+nkpt+9*nsym)-8*15 !non allocatable var. 2084 ! 2085 write(msg,'(2a,f8.2,a)')ch10,& 2086 & ' Additional amount of memory required by "outkss" routine=',memsize*b2Mb,' Mbytes.' 2087 call wrtout(std_out,msg,'COLL') 2088 ! 2089 end subroutine memkss
m_io_kss/k2gamma_centered [ Functions ]
[ Top ] [ m_io_kss ] [ Functions ]
NAME
k2gamma_centered
FUNCTION
Helper function to translate a set of wavefunctions from the k-centered G-sphere to the Gamma-centered G-sphere used for GW calculations.
INPUTS
npw_k=Number of planewaves in the k-centered basis set. kss_npw=Number of planewaves in the Gamma-centered G-sphere. nspinor=Number of spinorial component. nbandksseff=Number of bands in input-output arrays. [icg]=Shift to be used when accessing the cg array. 0 if not specified (usually k_index). [eig_vec(2,npw_k*nspinor,nbandksseff)]=wavefunctions defined on the k-centered G-sphere. [cg(2,ikg+1:ikg+npw_k*nspinor*nbandksseff)]=wavefunctions defined on the k-centered G-sphere. ngfft(18)=Info on the FFT. MPI_enreg<MPI_type>=Structure gathering info on the parallelization. istwf_k ecut gbig(3,kss_npw) kg_k(3,npw_k) gmet(3,3) kpoint(3)
OUTPUT
wfg(2,kss_npw*nspinor,nbandksseff)=Wavefunctions in the Gamma-centered representation.
NOTES
1) icg is used only if cg is present. 2) cg and eig_vec are mutually exclusive. One and only one can be passed to the routine.
SOURCE
595 subroutine k2gamma_centered(kpoint,npw_k,istwf_k,ecut,kg_k,kss_npw,nspinor,nbandksseff,ngfft,gmet,& 596 & MPI_enreg,gbig,ug,icg,cg,eig_vec) 597 598 !Arguments ------------------------------------ 599 !scalars 600 integer,intent(in) :: nbandksseff,nspinor,kss_npw,npw_k,istwf_k 601 integer,optional,intent(in) :: icg 602 real(dp),intent(in) :: ecut 603 type(MPI_type),intent(inout) :: MPI_enreg 604 !arrays 605 integer,intent(in) :: gbig(3,kss_npw) 606 integer,intent(in) :: kg_k(3,npw_k) 607 integer,intent(in) :: ngfft(18) 608 real(dp),intent(in) :: gmet(3,3),kpoint(3) 609 real(dp),intent(out) :: ug(2,kss_npw*nspinor,nbandksseff) 610 real(dp),optional,intent(in) :: eig_vec(2,npw_k*nspinor,nbandksseff) 611 real(dp),optional,intent(in) :: cg(:,:) 612 613 !Local variables------------------------------- 614 !scalars 615 integer,parameter :: tobox=1,tosph=-1 616 integer :: band,ispinor,spinor_shift2,spinor_shift1,ig,my_icg,ierr 617 integer :: n1,n2,n3,n4,n5,n6,ndat,full_npw_k,ii 618 character(len=500) :: msg 619 !arrays 620 integer :: identity(3,3)=RESHAPE((/1,0,0,0,1,0,0,0,1/),(/3,3/)) 621 integer :: no_shift(3)=(/0,0,0/) 622 integer,allocatable :: trsl(:),full_kg_k(:,:) 623 real(dp),allocatable :: cfft(:,:,:,:) 624 real(dp),allocatable :: full_cg(:,:),tmp_cg(:,:) 625 626 ! ********************************************************************* 627 628 if (PRESENT(cg).and.PRESENT(eig_vec)) then 629 ABI_ERROR("Both cg and eig_vec are present!") 630 end if 631 632 ! Mapping between the gamma-centered basis set and the k-centered one. 633 ! trsl(ig)=npw_k+1 if vector ig is not inside the k-centered G-sphere. 634 ABI_MALLOC(trsl,(kss_npw)) 635 636 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 637 n4=ngfft(4); n5=ngfft(5); n6=ngfft(6) 638 639 if (istwf_k==1) then ! Full k-centered G-sphere. 640 call table_gbig2kg(npw_k,kg_k,kss_npw,gbig,trsl,ierr) 641 if (ierr/=0.and.(kss_npw>=npw_k)) then 642 ABI_ERROR(' The set of G vectors is inconsistent') 643 end if 644 645 else ! Calculate full kg with istwf_k=1 then do the mapping. 646 call get_kg(kpoint,1,ecut,gmet,full_npw_k,full_kg_k) 647 648 call table_gbig2kg(full_npw_k,full_kg_k,kss_npw,gbig,trsl,ierr) 649 if (ierr/=0.and.(kss_npw>=npw_k)) then 650 ABI_ERROR(' The set of G vectors is inconsistent') 651 end if 652 end if 653 ! 654 ! Branching, depending on optional arguments. 655 if (PRESENT(cg)) then 656 my_icg=0; if (PRESENT(icg)) my_icg=icg 657 658 SELECT CASE (istwf_k) 659 660 CASE (1) 661 do band=1,nbandksseff 662 do ispinor=1,nspinor 663 spinor_shift1=(ispinor-1)*kss_npw 664 spinor_shift2=(ispinor-1)*npw_k 665 do ig=1,kss_npw ! Retrieve the correct components 666 if (trsl(ig)<=npw_k) then 667 ug(:,ig+spinor_shift1,band)=cg(:,trsl(ig)+spinor_shift2+(band-1)*npw_k*nspinor+my_icg) 668 else 669 ug(:,ig+spinor_shift1,band)=zero 670 end if 671 end do 672 end do 673 end do 674 675 CASE (2:9) 676 677 ABI_CHECK(nspinor==1,"nspinor/=1!") 678 ! 679 ! Convert input wfs from reduced to full G-sphere. 680 ndat=1 681 ABI_MALLOC(cfft,(2,n4,n5,n6*ndat)) 682 ABI_MALLOC(full_cg,(2,full_npw_k*ndat)) 683 ABI_MALLOC(tmp_cg,(2,npw_k*ndat)) 684 685 !write(std_out,*)"npw_k, full_kg_k",npw_k,full_npw_k 686 687 do band=1,nbandksseff 688 ii = (band-1)*npw_k 689 tmp_cg = cg(:,my_icg+ii+1:my_icg+ii+npw_k) 690 !write(776,*)"band= ",band,tmp_cg !cg(1:,my_icg+1+ii:my_icg+ii+npw_k) 691 692 call sphere(tmp_cg,ndat,npw_k,cfft,n1,n2,n3,n4,n5,n6,kg_k,istwf_k,tobox,MPI_enreg%me_g0,no_shift,identity,one) 693 694 call sphere(full_cg,ndat,full_npw_k,cfft,n1,n2,n3,n4,n5,n6,full_kg_k,1,tosph,MPI_enreg%me_g0,no_shift,identity,one) 695 !write(777,*)"band= ",band,full_cg(:,:) 696 697 do ig=1,kss_npw ! Retrieve the correct components 698 if (trsl(ig)<=full_npw_k) then 699 ug(:,ig,band)=full_cg(:,trsl(ig)) 700 else 701 ug(:,ig,band)=zero 702 end if 703 end do 704 end do !band 705 706 ABI_FREE(cfft) 707 ABI_FREE(tmp_cg) 708 ABI_FREE(full_cg) 709 710 CASE DEFAULT 711 ABI_BUG("Wrong istwf_k") 712 END SELECT 713 714 else if (PRESENT(eig_vec)) then 715 716 SELECT CASE (istwf_k) 717 718 CASE (1) 719 do band=1,nbandksseff 720 do ispinor=1,nspinor 721 spinor_shift1=(ispinor-1)*kss_npw 722 spinor_shift2=(ispinor-1)*npw_k 723 do ig=1,kss_npw ! Retrieve the correct components 724 if (trsl(ig)<=npw_k) then 725 ug(:,ig+spinor_shift1,band)=eig_vec(:,trsl(ig)+spinor_shift2,band) 726 else 727 ug(:,ig+spinor_shift1,band)=zero 728 end if 729 end do 730 end do 731 end do 732 733 CASE DEFAULT 734 write(msg,'(a,i0)')" Unsupported value for istwf_k: ",istwf_k 735 ABI_ERROR(msg) 736 END SELECT 737 738 else 739 ABI_ERROR("neither cg not eig_vec are in input") 740 end if 741 742 ABI_FREE(trsl) 743 if (allocated(full_kg_k)) then 744 ABI_FREE(full_kg_k) 745 end if 746 747 end subroutine k2gamma_centered
m_io_kss/kss_calc_vkb [ Functions ]
[ Top ] [ m_io_kss ] [ Functions ]
NAME
kss_calc_vkb
FUNCTION
This routine calculates the Kleynman-Bylander form factors and its derivatives needed for the evaluation of the matrix elements of the dipole operator <phi1|r|phi2>.
INPUTS
npw_k=Number of plane waves for this k-point. Psps<pseudopotential_type>=Structured datatype gathering information on the pseudopotentials. kg_k(3,npw_k)=Reduced coordinates of the G-vectors. kpoint(3)=The k-point in reduced coordinates. rprimd(3,3)=dimensional primitive translations for real space (bohr)
OUTPUT
vkb (npw_k,Psps%ntypat,Psps%mpsang)=KB form factors. vkbd(npw_k,Psps%ntypat,Psps%mpsang)=KB form factor derivatives. vkbsign(Psps%mpsang,Psps%ntypat) =KS dyadic sign.
NOTES
This piece of code has been extracted from outkss.F90. The implementation is consistent with the KSS file formata (Fortran version) but it presents two design flaws. 1) Pseudo with more that one projector per l-channel are not supported. 2) Ordering of dimensions in vkb and vkbd is not optimal. We are not programming C!!!
TODO
*) Spinorial case is not implemented.
SOURCE
934 subroutine kss_calc_vkb(Psps,kpoint,npw_k,kg_k,rprimd,vkbsign,vkb,vkbd) 935 936 !Arguments ------------------------------------ 937 !scalars 938 integer,intent(in) :: npw_k 939 type(Pseudopotential_type),intent(in) :: Psps 940 !arrays 941 integer,intent(in) :: kg_k(3,npw_k) 942 real(dp),intent(in) :: kpoint(3),rprimd(3,3) 943 real(dp),intent(out) :: vkb (npw_k,Psps%ntypat,Psps%mpsang) 944 real(dp),intent(out) :: vkbd(npw_k,Psps%ntypat,Psps%mpsang) 945 real(dp),intent(out) :: vkbsign(Psps%mpsang,Psps%ntypat) 946 947 !Local variables ------------------------------ 948 !scalars 949 integer :: dimffnl,ider,idir,itypat,nkpg,il0,in 950 integer :: il,ilmn,ig,is 951 real(dp) :: ucvol,effmass_free,ecutsm,ecut 952 !arrays 953 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 954 real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:) 955 real(dp),allocatable :: ylm(:,:),ylm_gr(:,:,:),ylm_k(:,:) 956 957 ! ************************************************************************* 958 959 DBG_ENTER("COLL") 960 961 ABI_CHECK(Psps%usepaw==0,"You should not be here!") 962 ABI_CHECK(Psps%useylm==0,"useylm/=0 not considered!") 963 964 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 965 ! 966 ! === Save KB dyadic sign (integer-valued) === 967 vkbsign=zero 968 do itypat=1,Psps%ntypat 969 il0=0 970 do ilmn=1,Psps%lmnmax 971 il=1+Psps%indlmn(1,ilmn,itypat) 972 in=Psps%indlmn(3,ilmn,itypat) 973 if (il/=il0 .and. in==1) then 974 il0=il 975 vkbsign(il,itypat)=DSIGN(one,Psps%ekb(ilmn,itypat)) 976 end if 977 end do 978 end do 979 980 ! === Allocate KB form factor and derivative wrt k+G === 981 ! * Here we do not use correct ordering for dimensions 982 983 ider=1; dimffnl=2 ! To retrieve the first derivative. 984 idir=0; nkpg=0 985 ! 986 ! Quantities used only if useylm==1 987 ABI_MALLOC(ylm,(npw_k,Psps%mpsang**2*Psps%useylm)) 988 ABI_MALLOC(ylm_gr,(npw_k,3+6*(ider/2),Psps%mpsang**2*Psps%useylm)) 989 ABI_MALLOC(ylm_k,(npw_k,Psps%mpsang**2*Psps%useylm)) 990 ABI_MALLOC(kpg_dum,(npw_k,nkpg)) 991 992 ABI_MALLOC(ffnl,(npw_k,dimffnl,Psps%lmnmax,Psps%ntypat)) 993 994 call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,gmet,gprimd,ider,idir,Psps%indlmn,& 995 kg_k,kpg_dum,kpoint,Psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,& 996 Psps%ntypat,Psps%pspso,Psps%qgrid_ff,rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr) 997 998 ABI_FREE(kpg_dum) 999 ABI_FREE(ylm) 1000 ABI_FREE(ylm_gr) 1001 ABI_FREE(ylm_k) 1002 1003 ABI_MALLOC(modkplusg,(npw_k)) 1004 1005 effmass_free=one; ecutsm=zero; ecut=HUGE(one) 1006 ! call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,modkplusg,kpoint,npw_k) 1007 call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,modkplusg,kpoint,npw_k,0,0) 1008 modkplusg(:)=SQRT(half/pi**2*modkplusg(:)) 1009 modkplusg(:)=MAX(modkplusg(:),tol10) 1010 1011 !do ig=1,npw_k 1012 ! kpg(:)= kpoint(:)+kg_k(:,ig) 1013 ! modkplusg(ig) = normv(kpg,gmet,"G") 1014 !end do 1015 1016 ! Calculate matrix elements. 1017 vkb=zero; vkbd=zero 1018 1019 do is=1,Psps%ntypat 1020 il0=0 1021 do ilmn=1,Psps%lmnmax 1022 il=1+Psps%indlmn(1,ilmn,is) 1023 in=Psps%indlmn(3,ilmn,is) 1024 if ((il/=il0).and.(in==1)) then 1025 il0=il 1026 if (ABS(Psps%ekb(ilmn,is))>1.0d-10) then 1027 if (il==1) then 1028 vkb (1:npw_k,is,il) = ffnl(:,1,ilmn,is) 1029 vkbd(1:npw_k,is,il) = ffnl(:,2,ilmn,is)*modkplusg(:)/two_pi 1030 else if (il==2) then 1031 vkb(1:npw_k,is,il) = ffnl(:,1,ilmn,is)*modkplusg(:) 1032 do ig=1,npw_k 1033 vkbd(ig,is,il) = ((ffnl(ig,2,ilmn,is)*modkplusg(ig)*modkplusg(ig))+& 1034 ffnl(ig,1,ilmn,is) )/two_pi 1035 end do 1036 else if (il==3) then 1037 vkb (1:npw_k,is,il) = ffnl(:,1,ilmn,is)*modkplusg(:)**2 1038 vkbd(1:npw_k,is,il) = (ffnl(:,2,ilmn,is)*modkplusg(:)**3+& 1039 2*ffnl(:,1,ilmn,is)*modkplusg(:) )/two_pi 1040 else if (il==4) then 1041 vkb (1:npw_k,is,il) = ffnl(:,1,ilmn,is)*modkplusg(:)**3 1042 vkbd(1:npw_k,is,il) = (ffnl(:,2,ilmn,is)*modkplusg(:)**4+& 1043 3*ffnl(:,1,ilmn,is)*modkplusg(:)**2 )/two_pi 1044 end if 1045 vkb (:,is,il) = SQRT(4*pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,is)))*vkb (:,is,il) 1046 vkbd(:,is,il) = SQRT(4*pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,is)))*vkbd(:,is,il) 1047 else 1048 vkb (:,is,il)=zero 1049 vkbd(:,is,il)=zero 1050 end if 1051 end if 1052 end do 1053 end do 1054 1055 ABI_FREE(ffnl) 1056 ABI_FREE(modkplusg) 1057 1058 DBG_EXIT("COLL") 1059 1060 end subroutine kss_calc_vkb
m_io_kss/make_gvec_kss [ Functions ]
[ Top ] [ m_io_kss ] [ Functions ]
NAME
make_gvec_kss
FUNCTION
Build the list of G-vectors using the KSS convention.
INPUTS
nkpt=Number of k-points.
nsym=Number of symmetries.
prtvol=Verbosity option.
symmorphi=
0 : Old (Obsolete) implementation => Suppress inversion from symmetries list
1 : Use input symrel, tnons.
ecut_eff=Effective cutoff
symrel(3,3,nsym)= Symmetry operation in real space.
tnons(3,nsym)=Fractional translations
kptns(3,nkpt)=K-points in reduced coordinates.
OUTPUT
npwkss = Input: Initial guess for the number of G-vectors required. Use 0 to have the
full list of G-vectors that form a closed shell.
Output: Actual number of G-vectors that form a set of closed shells
gvec_kss(:,:) = Input: null pointer. Output: gvec_kss(3,npwkss), list of G-vectors (closed shells)
ierr=Status error
SOURCE
780 subroutine make_gvec_kss(nkpt,kptns,ecut_eff,symmorphi,nsym,symrel,tnons,gprimd,prtvol,npwkss,gvec_kss,ierr) 781 782 !Arguments ------------------------------------ 783 !scalars 784 integer,intent(in) :: nkpt,nsym,prtvol,symmorphi 785 integer,intent(out) :: ierr 786 integer,intent(inout) :: npwkss 787 real(dp),intent(in) :: ecut_eff 788 !arrays 789 integer,intent(in) :: symrel(3,3,nsym) 790 integer,pointer :: gvec_kss(:,:) 791 real(dp),intent(in) :: tnons(3,nsym),kptns(3,nkpt) 792 real(dp),intent(in) :: gprimd(3,3) 793 794 !Local variables------------------------------- 795 !scalars 796 integer :: ii,ishm,maxpw,nbase 797 integer :: nrst1,nrst2,nsym2,pinv 798 integer,pointer :: gbig(:,:) 799 character(len=500) :: msg 800 !arrays 801 integer,pointer :: symrel2(:,:,:),shlim(:) 802 real(dp),pointer :: tnons2(:,:) 803 ! ********************************************************************* 804 805 ierr = 0 806 write(msg,'(2a)')ch10,' Sorting g-vecs for an output of states on an unique "big" PW basis.' 807 call wrtout(std_out,msg,'COLL') 808 809 !ecut_eff = ecut * Dtset%dilatmx**2 ! Use ecut_eff instead of ecut_eff since otherwise 810 ! 811 !============================================================ 812 !=== Prepare set containing all G-vectors sorted by stars === 813 !============================================================ 814 ! 815 !=== Analyze symmetry operations === 816 if (symmorphi==0) then ! Old (Obsolete) implementation: Suppress inversion from symmetries list: 817 nullify(symrel2,tnons2) 818 call remove_inversion(nsym,symrel,tnons,nsym2,symrel2,tnons2,pinv) 819 if (ANY(ABS(tnons2(:,1:nsym2))>tol8)) then 820 write(msg,'(3a)')& 821 & ' Non-symmorphic operations still remain in the symmetries list ',ch10,& 822 & ' Program does not stop but _KSS file will not be created...' 823 ABI_WARNING(msg) 824 ierr=ierr+1 ; RETURN 825 end if 826 else if (symmorphi==1) then 827 ! If in the input file symmorphi==1 all the symmetry operations are retained: 828 ! both identity and inversion (if any) as well as non-symmorphic operations. 829 nsym2=nsym ; pinv=1 830 ABI_MALLOC(symrel2,(3,3,nsym)) 831 ABI_MALLOC(tnons2,(3,nsym)) 832 symrel2(:,:,:)=symrel(:,:,1:nsym) 833 tnons2(:,:) =tnons(:,1:nsym) 834 else 835 write(msg,'(a,i4,3a)')& 836 & ' symmorphi = ',symmorphi,' while it must be 0 or 1',ch10,& 837 & ' Program does not stop but KSS file will not be created...' 838 ABI_WARNING(msg) 839 ierr=ierr+1 ; RETURN 840 end if 841 ! 842 !=================================================================== 843 !==== Merge the set of k-centered G-spheres into a big set gbig ==== 844 !=================================================================== 845 !* Vectors in gbig are ordered by shells 846 ! 847 nullify(gbig,shlim) 848 call merge_and_sort_kg(nkpt,kptns,ecut_eff,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p=shlim) 849 850 nbase = SIZE(shlim) ! Number of independent G in the big sphere. 851 maxpw = shlim(nbase) ! Total number of G"s in the big sphere. 852 ! 853 ! * Determine optimal number of bands and G"s to be written. 854 !npwkss=Dtset%npwkss 855 if ((npwkss==0).or.(npwkss>=maxpw)) then 856 npwkss=maxpw 857 write(msg,'(5a)')& 858 & ' Since the number of g''s to be written on file',ch10,& 859 & ' was 0 or too large, it has been set to the max. value.,',ch10,& 860 & ' computed from the union of the sets of G vectors for the different k-points.' 861 call wrtout(std_out,msg,'COLL') 862 end if 863 864 ishm=0 865 do ii=1,nbase 866 if (shlim(ii)<=npwkss) then 867 ishm=ii 868 else 869 EXIT 870 end if 871 end do 872 !ishm=bisect(shlim,npwkss) 873 874 if (shlim(ishm)/=npwkss) then 875 nrst1=shlim(ishm) 876 nrst2=MIN0(shlim(MIN0(ishm+1,nbase)),maxpw) 877 if (IABS(npwkss-nrst2)<IABS(npwkss-nrst1)) nrst1=nrst2 878 npwkss=nrst1 879 if (shlim(ishm)<npwkss) ishm=ishm+1 880 write(msg,'(3a)')& 881 & ' The number of G''s to be written on file is not a whole number of stars ',ch10,& 882 & ' the program set it to the nearest star limit.' 883 call wrtout(std_out,msg,'COLL') 884 end if 885 886 write(msg,'(a,i5)')' Number of G-vectors is: ',npwkss 887 call wrtout(std_out,msg,'COLL') 888 889 ABI_MALLOC(gvec_kss,(3,npwkss)) 890 gvec_kss = gbig(:,1:npwkss) 891 892 ABI_FREE(gbig) 893 ABI_FREE(symrel2) 894 ABI_FREE(tnons2) 895 ABI_FREE(shlim) 896 897 end subroutine make_gvec_kss
m_io_kss/outkss [ Functions ]
[ Top ] [ m_io_kss ] [ Functions ]
NAME
outkss
FUNCTION
This routine creates an output file containing the Kohn-Sham electronic Structure
for a large number of eigenstates (energies and eigen-functions).
The resulting file (_KSS) is needed for a GW post-treatment.
The routine drives the following operations:
- Re-ordering G-vectors according to stars (sets of Gs related by symmetry operations).
A set of g for all k-points is created.
- Creating and opening the output "_KSS'" file
- Printing out output file header information...
... and, for each k-point:
According to 'kssform', either
- Re-computing <G|H|G_prim> matrix elements for all (G, G_prim).
Diagonalizing H in the plane-wave basis.
or - Taking eigenvalues/vectors from congugate-gradient ones.
- Writing out eigenvalues and eigenvectors.
INPUTS
cg(2,mcg)=planewave coefficients of wavefunctions. usecprj=1 if cprj datastructure has been allocated (ONLY PAW) Cprj(natom,mcprj*usecprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with all NL projectors (only for PAW) NOTE that Cprj are unsorted, see ctoprj.F90 Dtfil <type(datafiles_type)>=variables related to files Dtset <type(dataset_type)>=all input variables for this dataset ecut=cut-off energy for plane wave basis sphere (Ha) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=dimensional reciprocal space primitive translations Hdr <type(hdr_type)>=the header of wf, den and pot files kssform=govern the Kohn-Sham Structure file format mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfft=maximum size of 1D FFTs mkmem =number of k points treated by this node. MPI_enreg=information about MPI parallelization mpsang= 1+maximum angular momentum for nonlocal pseudopotentials mpw=maximum dimensioned size of npw. my_natom=number of atoms treated by current processor natom=number of atoms in cell. nfft=(effective) number of FFT grid points (for this processor) nkpt=number of k points. npwarr(nkpt)=number of planewaves in basis at this k point nsppol=1 for unpolarized, 2 for spin-polarized nspden=number of density components nsym=number of symmetries in space group ntypat=number of types of atoms in unit cell. occ(mband*nkpt*nsppol)=occupation number for each band (usually 2) for each k. Pawtab(Psps%ntypat*Psps%usepaw) <type(pawtab_type)>=paw tabulated starting data Pawfgr<pawfgr_type>=fine grid parameters and related data prtvol=control print volume and debugging output Psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translations for real space (bohr) vtrial(nfft,nspden)=the trial potential xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
Output is written on file. ierr=Status error.
NOTES
* This routine is maintained for legacy reasons. Abinit8 is not able to read KSS files
anymore hence KSS files should be used only to interface Abinit with external codes
that are still using the old KSS format.
* The routine can be time consuming (in particular when computing
<G|H|G_prim> elements for all (G, G_prim)) (kssform=1).
So, it is recommended to call it once per run...
* The IO code is not parallelized and this represents a serious bottleneck when np is large.
* when kssform==1, the routine RE-computes all Hamiltonian terms.
So it is equivalent to an additional electronic SC cycle.
(This has no effect is convergence was reach...
If not, eigenvalues/vectors may differs from the congugaste gradient ones)
* The KB form factors and derivatives are not calculated correctly if there are
pseudos with more than one projector in an angular momentum channel.
* In the ETSF output format (Dtset%iomode == 3), the complete symmetry set
is output. So, if reading programs need only the symmorphic symmetries, they
will need to remove themselves the non-symmorphic ones.
* There exists two file formats:
kssform==1 diagonalized file _KSS in real(dp) is generated.
kssform==3 same as kssform=1 but the wavefunctions are not diagonalized
(they are taken from conjugate-gradient ones)
Old kssform=0 and kssform=2 are obsolete and no longer available
TESTS
* ETSF_IO output is tested in tests/etsf_io/t02.
SOURCE
1163 subroutine outkss(crystal,Dtfil,Dtset,ecut,gmet,gprimd,Hdr,& 1164 & kssform,mband,mcg,mcprj,mgfft,mkmem,MPI_enreg,mpsang,mpw,my_natom,natom,& 1165 & nfft,nkpt,npwarr,nspden,nsppol,nsym,ntypat,occ,Pawtab,Pawfgr,Paw_ij,& 1166 & prtvol,Psps,rprimd,vtrial,xred,cg,usecprj,Cprj,eigen,ierr) 1167 1168 use m_linalg_interfaces 1169 1170 !Arguments ------------------------------------ 1171 !scalars 1172 integer,intent(in) :: kssform,mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,usecprj 1173 integer,intent(in) :: nfft,nkpt,nsppol,nspden,nsym,ntypat,prtvol 1174 integer,intent(out) :: ierr 1175 real(dp),intent(in) :: ecut 1176 type(MPI_type),intent(inout) :: MPI_enreg 1177 type(Datafiles_type),intent(in) :: Dtfil 1178 type(Dataset_type),intent(in) :: Dtset 1179 type(Hdr_type),intent(inout) :: Hdr 1180 type(Pseudopotential_type),intent(in) :: Psps 1181 type(pawfgr_type), intent(in) :: Pawfgr 1182 type(crystal_t),intent(in) :: crystal 1183 !arrays 1184 integer,intent(in),target :: npwarr(nkpt) 1185 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(mband*nkpt*nsppol) 1186 real(dp),intent(in) :: rprimd(3,3) 1187 real(dp),intent(inout) :: vtrial(nfft,nspden) 1188 real(dp),intent(in) :: xred(3,natom) 1189 real(dp),intent(in) :: cg(2,mcg),eigen(mband*nkpt*nsppol) 1190 type(pawcprj_type),intent(in) :: Cprj(natom,mcprj*usecprj) 1191 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 1192 type(paw_ij_type),intent(inout),target :: Paw_ij(my_natom*Psps%usepaw) 1193 1194 !Local variables------------------------------- 1195 !scalars 1196 integer,parameter :: tim_rwwf=0,bufnb=20 1197 integer :: untkss,onband_diago 1198 integer :: bdtot_index,i,ib,ibp,iomode 1199 integer :: ibsp,ibsp1,ibsp2,ibg,ig,ii,ikpt 1200 integer :: master,receiver,sender,spinor_shift1,shift 1201 integer :: ishm,ispinor,isppol,istwf_k,my_rank,j 1202 integer :: k_index,maxpw,mproj,mtag,n1,n2,n2dim,n3,n4,n5,n6,nband_k 1203 integer :: nbandkss_k,nbandksseff,nbase,nprocs,npw_k,onpw_k,npwkss 1204 integer :: nrst1,nrst2,nsym2,ntemp,pinv,sizepw,spaceComm,comm_self 1205 integer :: pad1,pad2 1206 integer :: bufrt,bufsz 1207 real(dp) :: cinf=1.0e24,csup=zero,einf=1.0e24,esup=zero 1208 real(dp) :: norm,cfact,ecut_eff 1209 logical :: do_diago,found,ltest,lhack 1210 logical,parameter :: skip_test_ortho=.FALSE. 1211 character(len=500) :: msg 1212 character(len=80) :: frmt1,frmt2 1213 character(len=10) :: stag(2)=(/' ',' '/) 1214 !arrays 1215 integer :: nbandkssk(nkpt) 1216 integer,pointer :: symrel2(:,:,:) 1217 integer,pointer :: gbig(:,:) 1218 integer,pointer :: shlim(:) 1219 integer,allocatable :: kg_k(:,:) 1220 integer,allocatable :: dimlmn(:) 1221 integer :: nattyp_dum(0) 1222 real(dp) :: ovlp(2),kpoint(3),tsec(2) 1223 real(dp),pointer :: tnons2(:,:) 1224 real(dp),allocatable :: ene(:) 1225 real(dp),pointer :: eig_ene(:),eig_vec(:,:,:) 1226 real(dp),allocatable :: occ_k(:) 1227 real(dp),allocatable,target :: wfg(:,:,:) 1228 real(dp),ABI_CONTIGUOUS pointer :: ug1(:,:),ug2(:,:) 1229 type(pawcprj_type),allocatable :: Cprjnk_k(:,:) 1230 type(pawcprj_type),pointer :: Cprj_diago_k(:,:) 1231 type(ddiago_ctl_type) :: Diago_ctl 1232 type(paw_ij_type),pointer :: Paw_ij_all(:) 1233 ! ********************************************************************* 1234 1235 ABI_UNUSED(mkmem) 1236 1237 DBG_ENTER("COLL") 1238 1239 call timab(933,1,tsec) ! outkss 1240 call timab(934,1,tsec) ! outkss(Gsort+hd) 1241 1242 spaceComm=MPI_enreg%comm_cell 1243 my_rank=xmpi_comm_rank(spaceComm) 1244 nprocs=xmpi_comm_size(spaceComm) 1245 master=0 1246 1247 iomode = Dtset%iomode 1248 nullify(eig_ene) 1249 nullify(eig_vec) 1250 nullify(Cprj_diago_k) 1251 1252 ! JB: Valgrind complains about non initialized value. Set to -1 so if an array 1253 ! should be allocated with this "unintialized value" it crashes 1254 onband_diago = -1 1255 1256 !MG: since in seq case MPI_enreg%proc_distrb is not defined 1257 !we hack a bit the data type in order to get rid of MPI preprocessing options. 1258 !The previous status of %proc_distrb is restored before exiting. 1259 !Note that in case of seq run MPI_enreg%proc_distrb is nullified at the very beginning of abinit.F90 1260 ! 1261 !FIXME this is a design flaw that should be solved: proc_distrb should always 1262 !be allocated and filled with my_rank in case of sequential run otherwise checks like 1263 !if (nprocs>1.and.MPI_enreg%proc_distrb(ii)==me) leads to SIGFAULT under gfortran. 1264 !as the second array is not allocated. 1265 lhack=.FALSE. 1266 if (nprocs==1) then 1267 ltest=allocated(MPI_enreg%proc_distrb) 1268 if (.not.ltest) then 1269 ABI_MALLOC(MPI_enreg%proc_distrb,(nkpt,mband,nsppol)) 1270 MPI_enreg%proc_distrb=my_rank 1271 lhack=.TRUE. 1272 end if 1273 ltest=ALL(MPI_enreg%proc_distrb==my_rank) 1274 ABI_CHECK(ltest,'wrong values in %proc_distrb') 1275 end if 1276 ! 1277 !============================ 1278 !==== Perform some tests ==== 1279 !============================ 1280 ierr=0 1281 1282 if (iomode==IO_MODE_ETSF) then 1283 write(msg,'(3a)')& 1284 & 'when iomode==3 in outkss, support for netcdf ',ch10,& 1285 & 'must be compiled. Use --enable-netcdf when configuring ' 1286 #ifndef HAVE_NETCDF 1287 ABI_WARNING(msg) 1288 ierr = ierr + 1 1289 #endif 1290 end if 1291 1292 if (kssform==3) then 1293 write(msg,'(a,70("="),4a)')ch10,ch10,& 1294 & ' Calculating and writing out Kohn-Sham electronic Structure file',ch10, & 1295 & ' Using conjugate gradient wavefunctions and energies (kssform=3)' 1296 else if (kssform==1) then 1297 write(msg,'(a,70("="),4a,i1,a)') ch10,ch10, & 1298 & ' Calculating and writing out Kohn-Sham electronic Structure file',ch10, & 1299 & ' Using diagonalized wavefunctions and energies (kssform=',kssform,')' 1300 else 1301 write(msg,'(a,i0,2a)')& 1302 & " Unsupported value for kssform: ",kssform,ch10,& 1303 & " Program does not stop but _KSS file will not be created..." 1304 ierr=ierr+1 1305 end if 1306 call wrtout(std_out,msg,'COLL') 1307 call wrtout(ab_out,msg,'COLL') 1308 ! 1309 !* Check whether nband is constant in metals 1310 if ( (Dtset%occopt>=2.and.Dtset%occopt<=8) .and. (ANY(Dtset%nband(1:nkpt*nsppol)/=Dtset%nband(1))) ) then 1311 write(msg,'(3a,i4,a,i3,a,i4,3a)')& 1312 & ' The number of bands must be the same for all k-points ',ch10,& 1313 & ' but nband(1)=',Dtset%nband(1),' is different of nband(',& 1314 & ikpt+(isppol-1)*nkpt,')=',Dtset%nband(ikpt+(isppol-1)*nkpt),'.',ch10,& 1315 & ' Program does not stop but _KSS file will not be created...' 1316 ABI_WARNING(msg) 1317 ierr=ierr+1 1318 end if 1319 !* istwfk must be 1 for each k-point 1320 if (ANY(Dtset%istwfk(1:nkpt)/=1).and.kssform/=3) then 1321 write(msg,'(7a)')& 1322 & ' istwfk/=1 not allowed when kssform/=3 :',ch10,& 1323 & ' States output not programmed for time-reversal symmetry.',ch10,& 1324 & ' Action : change istwfk in input file (put it to 1 for all kpt).',ch10,& 1325 & ' Program does not stop but _KSS file will not be created...' 1326 ABI_WARNING(msg) 1327 ierr=ierr+1 1328 end if 1329 !* Check spin-orbit 1330 if (Psps%mpssoang/=mpsang) then 1331 write(msg,'(3a)')& 1332 & ' Variable mpspso should be 1 !',ch10,& 1333 & ' Program does not stop but _KSS file will not be created...' 1334 ABI_WARNING(msg) 1335 ierr=ierr+1 1336 end if 1337 !* Check mproj 1338 mproj=MAXVAL(Psps%indlmn(3,:,:)) 1339 if (mproj>1.and.Psps%usepaw==0) then ! TODO One has to derive the expression for [Vnl,r], in particular HGH and GTH psps 1340 write(msg,'(8a)')ch10,& 1341 & ' outkss : COMMENT - ',ch10,& 1342 & ' At least one NC pseudopotential has more that one projector per angular channel',ch10,& 1343 & ' Note that inclvkb==0 should be used in screening, since the evaluation of the commutator',ch10,& 1344 & ' for this particular case is not implemented yet' 1345 call wrtout(std_out,msg,'COLL') 1346 call wrtout(ab_out,msg,'COLL') 1347 end if 1348 !* Check max angular momentum 1349 if (MAXVAL(Psps%indlmn(1,:,:))+1 >= 5) then 1350 write(msg,'(3a)')& 1351 & ' Pseudopotentials with f-projectors not implemented',ch10,& 1352 & ' Program does not stop but _KSS file will not be created...' 1353 ABI_WARNING(msg) 1354 ierr=ierr+1 1355 end if 1356 !* Check useylm 1357 if (Psps%useylm/=0.and.Psps%usepaw==0) then 1358 write(msg,'(3a)')& 1359 & ' The present version of outkss does not work with useylm/=0 !',ch10,& 1360 & ' Program does not stop but _KSS file will not be created...' 1361 ABI_WARNING(msg) 1362 ierr=ierr+1 1363 end if 1364 !* Check PAW and kssform value 1365 if (Psps%usepaw/=0) then 1366 if (nprocs>1.and.kssform==1) then 1367 write(msg,'(3a)')& 1368 & ' Parallel PAW with kssform=1, not yet allowed',ch10,& 1369 & ' Program does not stop but _KSS file will not be created...' 1370 ABI_WARNING(msg) 1371 ierr=ierr+1 1372 end if 1373 if (kssform==3.and.usecprj/=1) then 1374 write(msg,'(3a)')& 1375 & ' If PAW and kssform=3, usecprj must be 1',ch10,& 1376 & ' Program does not stop but _KSS file will not be created...' 1377 ABI_WARNING(msg) 1378 ierr=ierr+1 1379 end if 1380 end if 1381 !* Check parallelization 1382 if (MPI_enreg%paralbd/=0) then 1383 write(msg,'(3a)')& 1384 & ' outkss cannot be used with parallelization on bands (paralbd/=0) !',ch10,& 1385 & ' Program does not stop but _KSS file will not be created...' 1386 ABI_WARNING(msg) 1387 ierr=ierr+1 1388 end if 1389 if (MPI_enreg%paral_spinor/=0) then 1390 write(msg,'(3a)')& 1391 & ' outkss cannot be used yet with parallelization on nspinors !',ch10,& 1392 & ' Program does not stop but _KSS file will not be created...' 1393 ABI_WARNING(msg) 1394 ierr=ierr+1 1395 1396 endif 1397 if (ierr/=0) then 1398 write(msg,'(3a)')& 1399 & ' outkss: Not allowed options found !',ch10,& 1400 & ' Program does not stop but _KSS file will not be created...' 1401 call wrtout(std_out,msg,'COLL') 1402 call wrtout(ab_out,msg,'COLL') 1403 write(msg,'(a)')' outkss: see the log file for more information.' 1404 call wrtout(ab_out,msg,'COLL') 1405 RETURN ! Houston we have a problem! 1406 end if 1407 ! 1408 !Estimate required memory in case of diagonalization. 1409 !TODO to be modified to take into account the case nsppol=2 1410 if (kssform/=3) then 1411 call memkss(mband,mgfft,mproj,Psps%mpssoang,mpw,natom,Dtset%ngfft,nkpt,dtset%nspinor,nsym,ntypat) 1412 end if 1413 ! 1414 !=== Initialize some variables === 1415 if (nsppol==2) stag(:)=(/'SPIN UP: ','SPIN DOWN:'/) 1416 n1=Dtset%ngfft(1); n2=Dtset%ngfft(2); n3=Dtset%ngfft(3) 1417 n4=Dtset%ngfft(4); n5=Dtset%ngfft(5); n6=Dtset%ngfft(6) 1418 ecut_eff = ecut * Dtset%dilatmx**2 ! Use ecut_eff instead of ecut_eff since otherwise 1419 !one cannot restart from a previous density file 1420 sizepw=2*mpw ; do_diago=(kssform/=3) 1421 ABI_MALLOC(dimlmn,(natom*Psps%usepaw)) 1422 if (Psps%usepaw==1) then 1423 call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,Dtset%typat,pawtab,'R') 1424 end if 1425 ! 1426 !============================================================ 1427 !=== Prepare set containing all G-vectors sorted by stars === 1428 !============================================================ 1429 write(msg,'(2a)')ch10,' Sorting g-vecs for an output of states on an unique "big" PW basis.' 1430 call wrtout(std_out,msg,'COLL') 1431 ! 1432 !=== Analyze symmetry operations === 1433 if (Dtset%symmorphi==0) then ! Old (Obsolete) implementation: Suppress inversion from symmetries list: 1434 nullify(symrel2,tnons2) 1435 call remove_inversion(nsym,Dtset%symrel,Dtset%tnons,nsym2,symrel2,tnons2,pinv) 1436 if (ANY(ABS(tnons2(:,1:nsym2))>tol8)) then 1437 write(msg,'(3a)')& 1438 & ' Non-symmorphic operations still remain in the symmetries list ',ch10,& 1439 & ' Program does not stop but _KSS file will not be created...' 1440 ABI_WARNING(msg) 1441 ierr=ierr+1 ; RETURN 1442 end if 1443 else if (Dtset%symmorphi==1) then 1444 ! If in the input file symmorphi==1 all the symmetry operations are retained: 1445 ! both identity and inversion (if any) as well as non-symmorphic operations. 1446 nsym2=nsym ; pinv=1 1447 ABI_MALLOC(symrel2,(3,3,nsym)) 1448 ABI_MALLOC(tnons2,(3,nsym)) 1449 symrel2(:,:,:)=Dtset%symrel(:,:,1:nsym) 1450 tnons2(:,:) =Dtset%tnons(:,1:nsym) 1451 else 1452 write(msg,'(a,i4,3a)')& 1453 & ' symmorphi = ',Dtset%symmorphi,' while it must be 0 or 1',ch10,& 1454 & ' Program does not stop but KSS file will not be created...' 1455 ABI_WARNING(msg) 1456 ierr=ierr+1 ; RETURN 1457 end if 1458 ! 1459 !=================================================================== 1460 !==== Merge the set of k-centered G-spheres into a big set gbig ==== 1461 !=================================================================== 1462 !* Vectors in gbig are ordered by shells 1463 ! 1464 nullify(gbig,shlim) 1465 call merge_and_sort_kg(nkpt,Dtset%kptns,ecut_eff,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p=shlim) 1466 1467 nbase = SIZE(shlim) ! Number of independent G in the big sphere. 1468 maxpw = shlim(nbase) ! Total number of G"s in the big sphere. 1469 ! 1470 !* Determine optimal number of bands and G"s to be written. 1471 npwkss=Dtset%npwkss 1472 if ((npwkss==0).or.(npwkss>=maxpw)) then 1473 npwkss=maxpw 1474 write(msg,'(5a)')& 1475 & ' Since the number of g''s to be written on file',ch10,& 1476 & ' was 0 or too large, it has been set to the max. value.,',ch10,& 1477 & ' computed from the union of the sets of G vectors for the different k-points.' 1478 call wrtout(std_out,msg,'COLL') 1479 end if 1480 1481 ishm=0 1482 do ii=1,nbase 1483 if (shlim(ii)<=npwkss) then 1484 ishm=ii 1485 else 1486 EXIT 1487 end if 1488 end do 1489 1490 if (shlim(ishm)/=npwkss) then 1491 nrst1=shlim(ishm) 1492 nrst2=MIN0(shlim(MIN0(ishm+1,nbase)),maxpw) 1493 if (IABS(npwkss-nrst2)<IABS(npwkss-nrst1)) nrst1=nrst2 1494 npwkss=nrst1 1495 if (shlim(ishm)<npwkss) ishm=ishm+1 1496 write(msg,'(3a)')& 1497 & ' The number of G''s to be written on file is not a whole number of stars ',ch10,& 1498 & ' the program set it to the nearest star limit.' 1499 call wrtout(std_out,msg,'COLL') 1500 end if 1501 1502 write(msg,'(a,i5)')' Number of g-vectors written on file is: ',npwkss 1503 call wrtout(std_out,msg,'COLL') 1504 ! 1505 !=== Check on the number of stored bands === 1506 if (do_diago) then 1507 1508 if (Dtset%nbandkss==-1.or.Dtset%nbandkss>=maxpw) then 1509 nbandkssk(1:nkpt)=npwarr(1:nkpt) 1510 write(msg,'(6a)')ch10,& 1511 & ' Since the number of bands to be computed was (-1) or',ch10,& 1512 & ' too large, it has been set to the max. value. allowed for each k,',ch10,& 1513 & ' thus, the minimum of the number of plane waves for each k point.' 1514 call wrtout(std_out,msg,'COLL') 1515 else 1516 nbandkssk(1:nkpt)=Dtset%nbandkss 1517 found=.FALSE. 1518 do ikpt=1,nkpt 1519 if (Dtset%nbandkss>npwarr(ikpt)) then 1520 nbandkssk(ikpt)=npwarr(ikpt) 1521 found=.TRUE. 1522 end if 1523 end do 1524 if (found) then 1525 write(msg,'(7a)')& 1526 & ' The value choosen for the number of bands in file',ch10,& 1527 & ' (nbandkss) was greater than at least one number of plane waves ',ch10,& 1528 & ' for a given k-point (npw_k).',ch10,' It has been modified consequently.' 1529 ABI_WARNING(msg) 1530 end if 1531 end if 1532 found=.FALSE. 1533 do ikpt=1,nkpt 1534 if (nbandkssk(ikpt)>npwkss) then 1535 nbandkssk(ikpt)=npwkss 1536 found=.TRUE. 1537 end if 1538 end do 1539 if (found) then 1540 write(msg,'(5a)')& 1541 & ' The number of bands to be computed (for one k) was',ch10,& 1542 & ' greater than the number of g-vectors to be written.',ch10,& 1543 & ' It has been modified consequently.' 1544 ABI_WARNING(msg) 1545 end if 1546 nbandksseff=MINVAL(nbandkssk) 1547 1548 else ! .not. do_diago 1549 do ikpt=1,nkpt 1550 do isppol=1,nsppol 1551 nbandkssk(ikpt)=Dtset%nband(ikpt+(isppol-1)*nkpt) 1552 end do 1553 end do 1554 nbandksseff=MINVAL(nbandkssk) 1555 if (Dtset%nbandkss>0 .and. Dtset%nbandkss<nbandksseff) then 1556 write(msg,'(a,i5,a,i5,2a)')& 1557 & ' Number of bands calculated=',nbandksseff,', greater than nbandkss=',Dtset%nbandkss,ch10,& 1558 & ' will write nbandkss bands on the KSS file' 1559 ABI_COMMENT(msg) 1560 nbandksseff=Dtset%nbandkss 1561 end if 1562 end if 1563 1564 write(msg,'(a,i5)')' Number of bands written on file is: ',nbandksseff 1565 call wrtout(std_out,msg,'COLL') 1566 1567 found= ANY(nbandkssk(1:nkpt)<npwarr(1:nkpt)) 1568 1569 if (do_diago) then 1570 if (found) then 1571 write(msg,'(6a)')ch10,& 1572 & ' Since the number of bands to be computed',ch10,& 1573 & ' is less than the number of G-vectors found,',ch10,& 1574 & ' the program will perform partial diagonalizations.' 1575 else 1576 write(msg,'(6a)')ch10,& 1577 & ' Since the number of bands to be computed',ch10,& 1578 & ' is equal to the nb of G-vectors found for each k-pt,',ch10,& 1579 & ' the program will perform complete diagonalizations.' 1580 end if 1581 call wrtout(std_out,msg,'COLL') 1582 end if 1583 ! 1584 !========================================================================== 1585 !=== Open KSS file for output, write header with dimensions and kb sign === 1586 !========================================================================== 1587 ! 1588 !* Output required disk space. 1589 call dsksta(ishm,Psps%usepaw,nbandksseff,mpsang,natom,ntypat,npwkss,nkpt,dtset%nspinor,nsppol,nsym2,dimlmn) 1590 1591 if (my_rank==master) then 1592 call write_kss_header(dtfil%fnameabo_kss,npwkss,ishm,nbandksseff,mband,nsym2,symrel2,tnons2,occ,gbig,shlim,& 1593 & crystal,Dtset,Hdr,Psps,iomode,untkss) 1594 end if 1595 1596 ABI_FREE(shlim) 1597 1598 if ( do_diago) msg = ' Diagonalized eigenvalues' 1599 if (.not.do_diago) msg = ' Conjugate gradient eigenvalues' 1600 call wrtout(ab_out,msg,'COLL') 1601 1602 if (Dtset%enunit==1) then 1603 msg=' k eigenvalues [eV]' 1604 else 1605 msg=' k eigenvalues [Hartree]' 1606 end if 1607 call wrtout(ab_out,msg,'COLL') 1608 ! 1609 !=== In case of PAW distributed atomic sites, need to retrieve the full paw_ij%dij === 1610 if (do_diago.and.Psps%usepaw==1.and.MPI_enreg%nproc_atom>1) then 1611 ABI_MALLOC(Paw_ij_all,(dtset%natom)) 1612 call paw_ij_gather(Paw_ij,Paw_ij_all,-1,MPI_enreg%comm_atom) 1613 else 1614 paw_ij_all => paw_ij 1615 end if 1616 1617 1618 call timab(934,2,tsec) ! outkss(Gsort+hd) 1619 ! 1620 1621 k_index=0; bdtot_index=0; ibg=0 1622 1623 do isppol=1,nsppol ! Loop over spins 1624 ! 1625 do ikpt=1,nkpt ! Loop over k-points. 1626 call timab(935,1,tsec) ! outkss(k-Loop) 1627 1628 nband_k =Dtset%nband(ikpt+(isppol-1)*nkpt) 1629 npw_k =npwarr(ikpt) 1630 istwf_k =Dtset%istwfk(ikpt) 1631 kpoint =Dtset%kptns(:,ikpt) 1632 nbandkss_k=nbandkssk(ikpt) 1633 mtag =5*(ikpt+(isppol-1)*nkpt) 1634 1635 1636 ! Get G-vectors, for this k-point. 1637 call get_kg(kpoint,istwf_k,ecut_eff,gmet,onpw_k,kg_k) 1638 ABI_CHECK(onpw_k==npw_k,"Mismatch in npw_k") 1639 ! 1640 ! ============================================ 1641 ! ==== Parallelism over k-points and spin ==== 1642 ! ============================================ 1643 if (MPI_enreg%proc_distrb(ikpt,1,isppol)==my_rank) then 1644 1645 write(msg,'(2a,i3,3x,a)')ch10,' k-point ',ikpt,stag(isppol) 1646 call wrtout(std_out, msg) 1647 1648 if (do_diago) then 1649 ! Direct diagonalization of the KS Hamiltonian. 1650 ABI_SFREE_PTR(eig_ene) 1651 ABI_SFREE_PTR(eig_vec) 1652 comm_self = xmpi_comm_self 1653 1654 call timab(936,1,tsec) 1655 1656 call init_ddiago_ctl(Diago_ctl,"Vectors",isppol,dtset%nspinor,ecut_eff,Dtset%kptns(:,ikpt),Dtset%nloalg,gmet,& 1657 nband_k=nbandkssk(ikpt),effmass_free=Dtset%effmass_free,istwf_k=Dtset%istwfk(ikpt),prtvol=Dtset%prtvol) 1658 1659 call ksdiago(Diago_ctl,nbandkssk(ikpt),Dtset%nfft,mgfft,Dtset%ngfft,natom,& 1660 Dtset%typat,nfft,dtset%nspinor,nspden,nsppol,Pawtab,Pawfgr,Paw_ij_all,& 1661 Psps,rprimd,vtrial,xred,onband_diago,eig_ene,eig_vec,Cprj_diago_k,comm_self,ierr) 1662 1663 call timab(936,2,tsec) 1664 end if 1665 1666 end if ! END of kpt+spin parallelism. 1667 ! 1668 ! =========================================================== 1669 ! ==== Transfer data between master and the working proc ==== 1670 ! =========================================================== 1671 call timab(937,1,tsec) !outkss(MPI_exch) 1672 ABI_MALLOC(Cprjnk_k,(0,0)) 1673 if (nprocs==1) then 1674 1675 if (Psps%usepaw==1) then ! Copy projectors for this k-point 1676 n2dim=min(nbandksseff*dtset%nspinor,onband_diago) 1677 if (kssform==3) n2dim=nband_k*dtset%nspinor 1678 ABI_FREE(Cprjnk_k) 1679 ABI_MALLOC(Cprjnk_k,(natom,n2dim)) 1680 call pawcprj_alloc(Cprjnk_k,0,dimlmn) 1681 if (kssform==3) then 1682 call pawcprj_copy(Cprj(:,ibg+1:ibg+dtset%nspinor*nband_k),Cprjnk_k) 1683 else 1684 !ABI_WARNING("Here I have to use onband_diago") !FIXME 1685 call pawcprj_copy(Cprj_diago_k(:,1:n2dim),Cprjnk_k) 1686 end if 1687 end if 1688 1689 else 1690 !parallel case 1691 1692 receiver=master; sender=MPI_enreg%proc_distrb(ikpt,1,isppol) 1693 1694 bufsz=nbandksseff/bufnb; bufrt=nbandksseff-bufnb*bufsz 1695 1696 if (my_rank==receiver.or.my_rank==sender) then 1697 1698 if (do_diago.and.(my_rank==receiver.and.my_rank/=sender)) then ! Alloc arrays if not done yet. 1699 ABI_MALLOC(eig_ene,(npw_k*dtset%nspinor)) 1700 ABI_MALLOC(eig_vec,(2,npw_k*dtset%nspinor,nbandkssk(ikpt))) 1701 end if 1702 1703 if (.not.do_diago) then 1704 1705 ABI_MALLOC(eig_vec,(2,npw_k*dtset%nspinor,nbandkssk(ikpt))) 1706 1707 if (my_rank==sender) then 1708 do ib=1,nbandksseff 1709 shift = k_index + (ib-1)*npw_k*dtset%nspinor 1710 do ig=1,npw_k*dtset%nspinor 1711 eig_vec(:,ig,ib)=cg(:,ig+shift) 1712 end do 1713 end do 1714 end if 1715 ! 1716 ! In case of PAW and kssform==3, retrieve matrix elements of the PAW projectors for this k-point 1717 if (Psps%usepaw==1) then 1718 n2dim=min(nbandksseff*dtset%nspinor,onband_diago) 1719 if (kssform==3) n2dim=nband_k*dtset%nspinor 1720 ABI_FREE(Cprjnk_k) 1721 ABI_MALLOC(Cprjnk_k,(natom,n2dim)) 1722 call pawcprj_alloc(Cprjnk_k,0,dimlmn) 1723 if (my_rank==sender) then 1724 if (kssform==3) then 1725 call pawcprj_copy(Cprj(:,ibg+1:ibg+dtset%nspinor*nband_k),Cprjnk_k) 1726 else 1727 !ABI_WARNING("Here I have to use onband_diago") !FIXME 1728 call pawcprj_copy(Cprj_diago_k(:,1:n2dim),Cprjnk_k) 1729 end if 1730 end if 1731 if (sender/=receiver) then 1732 call pawcprj_mpi_exch(natom,n2dim,dimlmn,0,Cprjnk_k,Cprjnk_k,sender,receiver,spaceComm,mtag+4,ierr) 1733 end if 1734 end if ! usepaw 1735 1736 else ! do_diago 1737 call xmpi_exch(eig_ene,nbandksseff,sender,eig_ene,receiver,spaceComm,mtag+1,ierr) 1738 end if 1739 1740 ! Exchange eigenvectors. 1741 if (bufsz>0) then 1742 do i=0,bufnb-1 1743 call xmpi_exch(eig_vec(:,:,i*bufsz+1:(i+1)*bufsz),2*npw_k*dtset%nspinor*bufsz,& 1744 & sender,eig_vec(:,:,i*bufsz+1:(i+1)*bufsz),receiver,spaceComm,mtag+2,ierr) 1745 end do 1746 end if 1747 if (bufrt>0) then 1748 call xmpi_exch(eig_vec(:,:,bufnb*bufsz+1:bufnb*bufsz+bufrt),2*npw_k*dtset%nspinor*bufrt,& 1749 & sender,eig_vec(:,:,bufnb*bufsz+1:bufnb*bufsz+bufrt),receiver,spaceComm,mtag+3,ierr) 1750 end if 1751 1752 end if 1753 end if !nprocs > 1 1754 call timab(937,2,tsec) !outkss(MPI_exch) 1755 1756 call timab(938,1,tsec) !outkss(write) 1757 1758 if (my_rank==master) then ! Prepare data for writing on disk. 1759 ABI_MALLOC(ene,(nbandksseff)) 1760 ABI_MALLOC(wfg,(2,npwkss*dtset%nspinor,nbandksseff)) 1761 ene=zero; wfg=zero 1762 1763 if (.not.do_diago) then 1764 ene(1:nbandksseff)=eigen(1+bdtot_index:nbandksseff+bdtot_index) 1765 1766 if (nprocs>1) then 1767 call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,& 1768 & MPI_enreg,gbig,wfg,eig_vec=eig_vec) 1769 else 1770 call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,& 1771 & MPI_enreg,gbig,wfg,icg=k_index,cg=cg) 1772 end if 1773 1774 else ! Direct diagonalization. 1775 ene(1:nbandksseff)=eig_ene(1:nbandksseff) 1776 1777 ! FIXME: receiver does not know Diago_ctl%npw_k 1778 call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,& 1779 & MPI_enreg,gbig,wfg,eig_vec=eig_vec) 1780 1781 ! * Check diagonalized eigenvalues with respect to conjugate gradient ones 1782 ntemp=MIN(nbandksseff,nband_k) 1783 if (ANY(ABS(ene(1:ntemp)-eigen(1+bdtot_index:ntemp+bdtot_index))>tol3)) then 1784 write(msg,'(3a)')& 1785 & ' The diagonalized eigenvalues differ by more than 10^-3 Hartree',ch10,& 1786 & ' with respect to the conjugated gradient values.' 1787 ABI_WARNING(msg) 1788 end if 1789 end if 1790 ! 1791 ! * Write out energies 1792 if (Dtset%enunit==1) then 1793 cfact=Ha_eV ; frmt1='(i4,4x,9(1x,f7.2))' ; frmt2='(8x,9(1x,f7.2))' 1794 write(msg,'(a,i3,3x,a)')' Eigenvalues in eV for ikpt= ',ikpt,stag(isppol) 1795 else 1796 cfact=one ; frmt1='(i4,4x,9(1x,f7.4))' ; frmt2='(8x,9(1x,f7.4))' 1797 write(msg,'(a,i3,3x,a)')' Eigenvalues in Hartree for ikpt= ',ikpt,stag(isppol) 1798 end if 1799 call wrtout(std_out,msg,'COLL') 1800 1801 write(msg,frmt1)ikpt,(ene(ib)*cfact,ib=1,MIN(9,nbandksseff)) 1802 call wrtout(std_out,msg,'COLL') 1803 call wrtout(ab_out,msg,'COLL') 1804 1805 if (nbandksseff>9) then 1806 do j=10,nbandksseff,9 1807 write(msg,frmt2) (ene(ib)*cfact,ib=j,MIN(j+8,nbandksseff)) 1808 call wrtout(std_out,msg,'COLL') 1809 call wrtout(ab_out,msg,'COLL') 1810 end do 1811 end if 1812 1813 if (skip_test_ortho) then ! Set this if to FALSE to skip test below 1814 einf=one; esup=one; cinf=zero; csup=zero 1815 else 1816 ! 1817 ! Test on the normalization of wavefunctions. 1818 ibsp=0 1819 do ib=1,nbandksseff 1820 norm=zero 1821 do ispinor=1,dtset%nspinor 1822 ibsp=ibsp+1 1823 spinor_shift1=(ispinor-1)*npwkss 1824 ug1 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ib) 1825 1826 !ovlp(1) =ddot(npwkss,ug1(1,:),1,ug1(1,:),1) + ddot(npwkss,ug1(2,:),1,ug1(2,:),1) 1827 ovlp(1) = cg_dznrm2(npwkss,ug1) 1828 ovlp(1) = ovlp(1)**2 1829 ovlp(2) = zero 1830 if (Psps%usepaw==1) ovlp = ovlp & 1831 & + paw_overlap(Cprjnk_k(:,ibsp:ibsp),Cprjnk_k(:,ibsp:ibsp),Dtset%typat,Pawtab,& 1832 & spinor_comm=MPI_enreg%comm_spinor) 1833 norm = norm + DABS(ovlp(1)) 1834 end do 1835 if (norm<einf) einf=norm 1836 if (norm>esup) esup=norm 1837 end do 1838 ! 1839 ! Test on the orthogonalization of wavefunctions. 1840 do ib=1,nbandksseff 1841 pad1=(ib-1)*dtset%nspinor 1842 do ibp=ib+1,nbandksseff 1843 pad2=(ibp-1)*dtset%nspinor 1844 ovlp(:)=zero 1845 do ispinor=1,dtset%nspinor 1846 ibsp1=pad1+ispinor 1847 ibsp2=pad2+ispinor 1848 spinor_shift1=(ispinor-1)*npwkss 1849 ug1 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ib ) 1850 ug2 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ibp) 1851 1852 !ovlp(1)=ddot(npwkss,ug1(1,:),1,ug2(1,:),1) + ddot(npwkss,ug1(2,:),1,ug2(2,:),1) 1853 !ovlp(2)=ddot(npwkss,ug1(1,:),1,ug2(2,:),1) - ddot(npwkss,ug1(2,:),1,ug2(1,:),1) 1854 ovlp = cg_zdotc(npwkss,ug1,ug2) 1855 1856 if (Psps%usepaw==1) ovlp= ovlp & 1857 & + paw_overlap(Cprjnk_k(:,ibsp1:ibsp1),Cprjnk_k(:,ibsp2:ibsp2),Dtset%typat,Pawtab,& 1858 & spinor_comm=MPI_enreg%comm_spinor) 1859 end do 1860 norm = DSQRT(ovlp(1)**2+ovlp(2)**2) 1861 if (norm<cinf) cinf=norm 1862 if (norm>csup) csup=norm 1863 end do 1864 end do 1865 end if 1866 1867 write(msg,'(a,i3,3x,a)')' Writing out eigenvalues/vectors for ikpt=',ikpt,stag(isppol) 1868 call wrtout(std_out,msg,'COLL') 1869 ! 1870 ! * Write occupation numbers on std_out. 1871 ABI_MALLOC(occ_k,(MAX(nband_k,nbandksseff))) 1872 occ_k(1:nband_k)=occ(1+bdtot_index:nband_k+bdtot_index) 1873 if (nband_k < nbandksseff) occ_k(nband_k+1:nbandksseff)=zero 1874 1875 write(msg,'(a,i3,3x,a)')' Occupation numbers for ikpt=',ikpt,stag(isppol) 1876 call wrtout(std_out,msg,'COLL') 1877 write(msg,'(i4,4x,9(1x,f7.4))')ikpt,(occ_k(ib),ib=1,MIN(9,nbandksseff)) 1878 call wrtout(std_out,msg,'COLL') 1879 if (nbandksseff>9) then 1880 do j=10,nbandksseff,9 1881 write(msg,'(8x,9(1x,f7.4))') (occ_k(ib),ib=j,MIN(j+8,nbandksseff)) 1882 call wrtout(std_out,msg,'COLL') 1883 end do 1884 end if 1885 ! 1886 ! ================================================================= 1887 ! ==== Write wavefunctions, KB and PAW matrix elements on disk ==== 1888 ! ================================================================= 1889 call write_kss_wfgk(untkss,ikpt,isppol,kpoint,dtset%nspinor,npwkss,& 1890 nbandksseff,natom,Psps,ene,occ_k,rprimd,gbig,wfg,Cprjnk_k,iomode) 1891 1892 ABI_FREE(occ_k) 1893 ABI_FREE(ene) 1894 ABI_FREE(wfg) 1895 1896 end if ! my_rank==master 1897 call timab(938,2,tsec) !outkss(write) 1898 1899 if (my_rank==master.or.my_rank==MPI_enreg%proc_distrb(ikpt,1,isppol)) then 1900 ABI_SFREE_PTR(eig_ene) 1901 ABI_SFREE_PTR(eig_vec) 1902 if (Psps%usepaw==1) call pawcprj_free(Cprjnk_k) 1903 end if 1904 ABI_FREE(Cprjnk_k) 1905 ABI_SFREE(kg_k) 1906 1907 ! if (MPI_enreg%paral_compil_kpt==1) then !cannot be used in seq run! 1908 if (.not.(proc_distrb_cycle(MPI_enreg%proc_distrb,ikpt,1,nband_k,isppol,my_rank))) then 1909 k_index=k_index+npw_k*nband_k*dtset%nspinor 1910 ibg=ibg+dtset%nspinor*nband_k 1911 end if 1912 bdtot_index=bdtot_index+nband_k 1913 1914 call xmpi_barrier(spaceComm) ! FIXME this barrier is detrimental in the case of direct diago! 1915 1916 call timab(935,2,tsec) !outkss(k-loop) 1917 end do ! ! End loop over k-points. 1918 end do ! spin 1919 1920 write(msg,'(3a,f9.6,2a,f9.6,4a,f9.6,2a,f9.6,a)')& 1921 & ' Test on the normalization of the wavefunctions',ch10,& 1922 & ' min sum_G |a(n,k,G)| = ',einf,ch10,& 1923 & ' max sum_G |a(n,k,G)| = ',esup,ch10,& 1924 & ' Test on the orthogonalization of the wavefunctions',ch10,& 1925 & ' min sum_G a(n,k,G)a(n'',k,G) = ',cinf,ch10,& 1926 & ' max sum_G a(n,k,G)a(n'',k,G) = ',csup,ch10 1927 call wrtout(std_out,msg,'COLL') 1928 call wrtout(ab_out,msg,'COLL') 1929 1930 ABI_FREE(gbig) 1931 ABI_FREE(symrel2) 1932 ABI_FREE(tnons2) 1933 if (Psps%usepaw==1) then 1934 ABI_FREE(dimlmn) 1935 if (do_diago.and.MPI_enreg%nproc_atom>1) then 1936 ABI_FREE(Paw_ij_all) 1937 end if 1938 end if 1939 ! 1940 !* Close file 1941 if (my_rank==master) then 1942 if (iomode==IO_MODE_FORTRAN) close(unit=untkss) 1943 #if defined HAVE_NETCDF 1944 if (iomode==IO_MODE_ETSF) then 1945 NCF_CHECK(nf90_close(untkss)) 1946 end if 1947 #endif 1948 end if 1949 1950 if (associated(Cprj_diago_k)) then 1951 call pawcprj_free(Cprj_diago_k) 1952 ABI_FREE(Cprj_diago_k) 1953 end if 1954 1955 if (lhack) then 1956 ABI_FREE(MPI_enreg%proc_distrb) 1957 end if 1958 1959 call wrtout(std_out, "outkss done", "COLL") 1960 call xmpi_barrier(spaceComm) 1961 1962 DBG_EXIT("COLL") 1963 call timab(933,2,tsec) ! outkss 1964 1965 contains
m_io_kss/write_kss_header [ Functions ]
[ Top ] [ m_io_kss ] [ Functions ]
NAME
write_kss_header
FUNCTION
Write the header of the KSS file either using plain Fortran-IO or netcdf with ETSF-IO format. Returns the unit number to be used for further writing. It should be executed by master node only.
INPUTS
filekss(len=fnlen)=The name of the KSS file. kss_npw=Number of planewaves used for the wavefunctions in the KSS files. ishm=Max number of shells written on file shlim(ishm)=The cumulative number of G"s in each shell. nbandksseff=Number of bands to be written. mband=The maximum number of bands treated by abinit. nsym2=Number of symmetry operations to be written on the header. symrel2(3,3,nsym2)=The symmetry operations in real space to be written. tnons2(3,nsym2)=The fractional translations associateed to symrel2. gbig(3,kss_npw)=The set of G-vectors for the KSS wavefunctions (Gamma-centered) Hdr<hdr_type>=The abinit header. Dtset <dataset_type>=all input variables for this dataset Psps<pseudopotential_type>=Structure gathering info on the pseudopotentials. iomode=Input variables specifying the fileformat. (0-->Fortran,3-->netcdf with ETSF-IO format). occ(mband*nkpt*nsppol)=The occupation factors.
OUTPUT
kss_unt=The unit number of the opened file.
SIDE EFFECTS
The KSS Header is written on file.
SOURCE
111 subroutine write_kss_header(filekss,kss_npw,ishm,nbandksseff,mband,nsym2,symrel2,tnons2,occ,gbig,shlim,& 112 crystal,Dtset,Hdr,Psps,iomode,kss_unt) 113 114 !Arguments ------------------------------------ 115 !scalars 116 integer,intent(in) :: iomode,kss_npw,nbandksseff,ishm,nsym2,mband 117 integer,intent(out) :: kss_unt 118 character(len=fnlen),intent(in) :: filekss 119 type(crystal_t),intent(in) :: crystal 120 type(pseudopotential_type),intent(in) :: Psps 121 type(Hdr_type),intent(in) :: Hdr 122 type(Dataset_type),intent(in) :: Dtset 123 !arrays 124 integer,intent(in) :: symrel2(3,3,nsym2) 125 integer,intent(in) :: gbig(3,kss_npw),shlim(ishm) 126 real(dp),intent(in) :: tnons2(3,nsym2) 127 real(dp),intent(in) :: occ(mband*Dtset%nkpt*Dtset%nsppol) 128 129 !Local variables------------------------------- 130 !scalars 131 integer :: nspinor,nsppol,nkpt,itypat,ierr 132 integer :: nb,isppol,ik,fform,ii,jj,kk,ig 133 integer :: il,il0,ilmn,in,ind1,ind2 134 character(len=80) :: title 135 character(len=500) :: msg 136 type(hdr_type) :: my_Hdr 137 type(dataset_type) :: Dtset_cpy 138 #ifdef HAVE_NETCDF 139 integer :: ncerr 140 #endif 141 !arrays 142 integer,allocatable :: vkbsign_int(:,:,:) 143 real(dp),allocatable :: vkbsign(:,:) 144 145 ! ********************************************************************* 146 147 DBG_ENTER("COLL") 148 149 nsppol = Dtset%nsppol 150 nkpt = Dtset%nkpt 151 nspinor= Dtset%nspinor 152 153 write(msg,'(3a)')ch10,' Opening file for KS structure output: ',TRIM(filekss) 154 call wrtout(std_out,msg,'COLL') 155 156 write(msg,'(a,i6)') ' number of Gamma centered plane waves ',kss_npw 157 call wrtout(std_out,msg,'COLL') 158 call wrtout(ab_out,msg,'COLL') 159 write(msg,'(a,i6)') ' number of Gamma centered shells ',ishm 160 call wrtout(std_out,msg,'COLL') 161 call wrtout(ab_out,msg,'COLL') 162 write(msg,'(a,i6)') ' number of bands ',nbandksseff 163 call wrtout(std_out,msg,'COLL') 164 call wrtout(ab_out,msg,'COLL') 165 write(msg,'(a,i6)') ' maximum angular momentum components ',Psps%mpsang 166 call wrtout(std_out,msg,'COLL') 167 call wrtout(ab_out,msg,'COLL') 168 write(msg,'(a,i2,a)')' number of symmetry operations ',nsym2,' (without inversion)' 169 call wrtout(std_out,msg,'COLL') 170 171 !Copy the header so that we can change some basic dimensions using the KSS values: 172 !(bantot, npwarr, nband) and the occupation factors 173 174 !Note that nsym and symrel might have been changed this has to be fixed 175 !carefully in the next patch since in the new implementation symmorphy=0 should be dafault 176 call hdr_copy(Hdr,my_Hdr) 177 178 my_Hdr%npwarr =kss_npw 179 my_Hdr%nband =nbandksseff 180 my_hdr%mband = maxval(my_hdr%nband) 181 my_Hdr%bantot =nbandksseff*nkpt*nsppol 182 183 my_Hdr%istwfk = 1 ! KSS file does not support istwfk/=1 even though the GS run 184 ! can take advantage of time-reversal symmetry. 185 186 !Copy the occ number in the new header with correct dimensions 187 !fill with zero the rest since mband can be < nbandksseff 188 !write(std_out,*)associated(my_Hdr%occ) 189 ABI_FREE(my_Hdr%occ) 190 ABI_MALLOC(my_Hdr%occ,(my_Hdr%bantot)) 191 !mband = MAXVAL(Hdr%nband) 192 193 my_Hdr%occ=zero; nb=MIN(mband,nbandksseff) 194 do isppol=1,nsppol 195 do ik=1,nkpt 196 ind1=1+(ik-1)*nbandksseff+(isppol-1)*nkpt*nbandksseff 197 ind2=1+(ik-1)*mband +(isppol-1)*nkpt*mband 198 my_Hdr%occ(ind1:ind1+nb-1) = occ(ind2:ind2+nb-1) 199 end do 200 end do 201 202 !Change dimension in the local Dtset_cpy as well. 203 dtset_cpy = Dtset%copy() 204 Dtset_cpy%mpw = kss_npw 205 Dtset_cpy%mband = nbandksseff 206 207 fform=502 208 209 SELECT CASE (iomode) 210 211 CASE (IO_MODE_FORTRAN) 212 213 if (open_file(filekss, msg, newunit=kss_unt, form="unformatted") /= 0) then 214 ABI_ERROR(msg) 215 end if 216 217 call my_hdr%fort_write(kss_unt, fform, ierr) 218 ABI_CHECK(ierr == 0, "hdr_Fort_write returned ierr != 0") 219 220 title='Results from ABINIT code'; write(kss_unt) title(1:80) 221 title='Ab-initio plane waves calculation'; write(kss_unt) title(1:80) 222 223 write(kss_unt) nsym2,nbandksseff,kss_npw,ishm,Psps%mpsang ! To be modified to deal with more than one projector 224 write(kss_unt) (((symrel2(ii,jj,kk),ii=1,3),jj=1,3),kk=1,nsym2) 225 write(kss_unt) ((tnons2(ii,kk),ii=1,3),kk=1,nsym2) 226 write(kss_unt) ((gbig(ii,ig),ii=1,3),ig=1,kss_npw) 227 write(kss_unt) (shlim(in),in=1,ishm) 228 229 ! Write vkbsign for NC pseudos with Fortran IO 230 ! MG FIXME: only one projector in each angular channel is treated. 231 ! Moreover the allocation is done in the wrong order for dimensions... 232 ! but if I change this code, compatibility with external codes is broken. 233 if (Psps%usepaw==0) then 234 ABI_MALLOC(vkbsign,(Psps%ntypat,Psps%mpsang)) 235 vkbsign(:,:)=zero 236 do itypat=1,Psps%ntypat 237 il0=0 238 do ilmn=1,Psps%lmnmax 239 il=1+Psps%indlmn(1,ilmn,itypat) 240 in=Psps%indlmn(3,ilmn,itypat) 241 if (il/=il0 .and. in==1) then 242 il0=il 243 vkbsign(itypat,il)=DSIGN(one,Psps%ekb(ilmn,itypat)) 244 end if 245 end do 246 end do 247 write(kss_unt) ((vkbsign(itypat,il),il=1,Psps%mpsang),itypat=1,Psps%ntypat) 248 ABI_FREE(vkbsign) 249 end if 250 251 #ifdef HAVE_NETCDF 252 CASE (IO_MODE_ETSF) 253 254 ! Create file. 255 NCF_CHECK(nctk_open_create(kss_unt, nctk_ncify(filekss), xmpi_comm_self)) 256 257 ! Add additional info from abinit header. 258 NCF_CHECK(my_hdr%ncwrite(kss_unt, fform, nc_define=.True.)) 259 260 ! Add info on crystalline structure 261 ! FIXME: Check symmorphi trick and crystal%symrel! 262 ! We currently use the dataset symmetries, as defined in the Hdr structure 263 ! instead of the symmetries recomputed in outkss. 264 NCF_CHECK(crystal%ncwrite(kss_unt)) 265 266 ! Defined G-vectors and wavefunctions. 267 call wfk_ncdef_dims_vars(kss_unt, my_hdr, fform, iskss=.True.) 268 !call abi_etsf_init(Dtset_cpy, filekss, 4, .false., my_Hdr%lmn_size, Psps, Dummy_wfs) 269 270 ! If NC pseudos, write vkbsign. 271 ! Here multi-projectors are supported, array is dimensioned according to etsf-io standard. 272 if (psps%usepaw == 0) then 273 274 ! Define dims and variables needed for KB matrix elements. 275 ncerr = nctk_def_dims(kss_unt, [ & 276 nctkdim_t("max_number_of_angular_momenta", psps%mpsang), & 277 nctkdim_t("max_number_of_projectors", psps%mproj) & 278 ]) 279 NCF_CHECK(ncerr) 280 281 ncerr = nctk_def_arrays(kss_unt, [ & 282 nctkarr_t("kb_formfactor_sign", "int", & 283 &"max_number_of_projectors, max_number_of_angular_momenta, number_of_atom_species"), & 284 nctkarr_t("kb_formfactors", "dp", & 285 &"max_number_of_coefficients, number_of_kpoints, max_number_of_projectors,& 286 &max_number_of_angular_momenta, number_of_atom_species"), & 287 nctkarr_t("kb_formfactor_derivative", "dp", & 288 &"max_number_of_coefficients, number_of_kpoints, max_number_of_projectors,& 289 &max_number_of_angular_momenta, number_of_atom_species") & 290 ]) 291 NCF_CHECK(ncerr) 292 293 ABI_MALLOC(vkbsign_int, (psps%mproj, Psps%mpsang, Psps%ntypat)) 294 vkbsign_int=0 295 do itypat=1,Psps%ntypat 296 do ilmn=1,Psps%lmnmax 297 il=1+Psps%indlmn(1,ilmn,itypat) 298 in=Psps%indlmn(3,ilmn,itypat) 299 vkbsign_int(in,il,itypat)=NINT(DSIGN(one,Psps%ekb(ilmn,itypat))) 300 end do 301 end do 302 303 NCF_CHECK(nctk_set_datamode(kss_unt)) 304 305 ! Write KB sign here 306 NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "kb_formfactor_sign"), vkbsign_int)) 307 ABI_FREE(vkbsign_int) 308 end if 309 310 NCF_CHECK(nctk_set_datamode(kss_unt)) 311 #endif 312 313 CASE DEFAULT 314 ABI_ERROR(sjoin("Unsupported value for iomode:", itoa(iomode))) 315 END SELECT 316 317 call Dtset_cpy%free() 318 call my_Hdr%free() 319 320 DBG_EXIT("COLL") 321 322 end subroutine write_kss_header
m_io_kss/write_kss_wfgk [ Functions ]
[ Top ] [ m_io_kss ] [ Functions ]
NAME
write_kss_wfgk
FUNCTION
Write the Gamma-centered wavefunctions and energies on the KSS file for a single k-point. (Only the master node should call this routine).
INPUTS
kss_unt=The unit number of the file ikpt=The index of the k-point isppol=The spin index. nspinor=number of spinorial components (on current proc) kss_npw=Number of planewaves used for the wavefunctions in the KSS files. npw_k=Number of plane-waves in the k-centered basis set. nbandksseff=Number of bands to be written. natom=Number of atoms. Psps<Pseudopotential_type>=Structure gathering pseudopotential data. kpoint(3)=The k-points in reduced coordinates. ene_k(nbandksseff)=Energies at this k-point occ_k(nbandksseff)=Occupation factors at this k-point. rprimd(3,3)=dimensional primitive translations for real space (bohr). gbig(3,kss_npw)=The set of G-vectors for the KSS wavefunctions (Gamma-centered) wfg(2,kss_npw*nspinor,nbandksseff)=The wavefunction Fourier coefficients. iomode=Input variables specifying the fileformat. (0-->Fortran,3--> netcdf with ETSF-IO format).
OUTPUT
Only writing.
SOURCE
472 subroutine write_kss_wfgk(kss_unt,ikpt,isppol,kpoint,nspinor,kss_npw,& 473 & nbandksseff,natom,Psps,ene_k,occ_k,rprimd,gbig,wfg,Cprjnk_k,iomode) 474 475 !Arguments ------------------------------------ 476 !scalars 477 integer,intent(in) :: ikpt,isppol,iomode,kss_npw,nspinor,kss_unt,nbandksseff 478 integer,intent(in) :: natom 479 type(pseudopotential_type),intent(in) :: Psps 480 !arrays 481 integer,intent(in) :: gbig(3,kss_npw) 482 real(dp),intent(in) :: kpoint(3),rprimd(3,3) 483 real(dp),intent(in) :: ene_k(nbandksseff),occ_k(nbandksseff) 484 real(dp),intent(in) :: wfg(2,kss_npw*nspinor,nbandksseff) 485 type(pawcprj_type),intent(in) :: Cprjnk_k(natom,nspinor*nbandksseff*Psps%usepaw) 486 487 !Local variables------------------------------- 488 !scalars 489 integer :: ib,ibsp,ig,ispinor,iatom,ii !,ierr 490 #ifdef HAVE_NETCDF 491 integer :: kg_varid,cg_varid,ncerr 492 character(len=nctk_slen) :: kdep 493 #endif 494 495 ! ********************************************************************* 496 497 ! Calculate and write KB form factors and derivative at this k-point. 498 if (Psps%usepaw==0) then 499 call write_vkb(kss_unt,ikpt,kpoint,kss_npw,gbig,rprimd,Psps,iomode) 500 end if 501 502 ! ============================================================ 503 ! ==== Write wavefunctions and PAW matrix elements on disk ==== 504 ! ============================================================ 505 SELECT CASE (iomode) 506 507 CASE (IO_MODE_FORTRAN) 508 write(kss_unt) (ene_k(ib),ib=1,nbandksseff) 509 510 ibsp=0 511 do ib=1,nbandksseff 512 write(kss_unt) (wfg(:,ig,ib),ig=1,kss_npw*nspinor) 513 if (Psps%usepaw==1) then ! Remember that cprj are unsorted. 514 do ispinor=1,nspinor 515 ibsp=ibsp+1 516 do iatom=1,natom 517 ii=Cprjnk_k(iatom,ibsp)%nlmn 518 write(kss_unt) (Cprjnk_k(iatom,ibsp)%cp(:,1:ii)) 519 end do 520 end do 521 end if 522 end do 523 524 #ifdef HAVE_NETCDF 525 CASE (IO_MODE_ETSF) 526 if (Psps%usepaw==1) then 527 ABI_WARNING("PAW output with ETSF-IO netcdf: cprj won't be written") 528 end if 529 530 ! Write G-vectors (gbig because it's not k-dependent) 531 NCF_CHECK(nf90_inq_varid(kss_unt, "reduced_coordinates_of_plane_waves", kg_varid)) 532 NCF_CHECK(nf90_get_att(kss_unt, kg_varid, "k_dependent", kdep)) 533 if (kdep == "no") then 534 ncerr = nf90_put_var(kss_unt, kg_varid, gbig, start=[1,1], count=[3,kss_npw]) 535 else 536 ncerr = nf90_put_var(kss_unt, kg_varid, gbig, start=[1,1,ikpt], count=[3,kss_npw,1]) 537 end if 538 NCF_CHECK_MSG(ncerr, "putting gibg") 539 540 ! Write wavefunctions 541 ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol] 542 NCF_CHECK(nf90_inq_varid(kss_unt, "coefficients_of_wavefunctions", cg_varid)) 543 ncerr = nf90_put_var(kss_unt, cg_varid, wfg, start=[1,1,1,1,ikpt,isppol], & 544 count=[2,kss_npw,nspinor,nbandksseff,1,1]) 545 NCF_CHECK_MSG(ncerr, "putting cg_k") 546 547 ! Write eigenvalues and occupations 548 NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "eigenvalues"), ene_k, start=[1,ikpt,isppol])) 549 NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "occupations"), occ_k, start=[1,ikpt,isppol])) 550 #endif 551 552 CASE DEFAULT 553 ABI_ERROR(sjoin("Unsupported iomode:", itoa(iomode))) 554 END SELECT 555 556 end subroutine write_kss_wfgk
m_io_kss/write_vkb [ Functions ]
[ Top ] [ m_io_kss ] [ Functions ]
NAME
write_vkb
FUNCTION
Writes the KB form factors and derivates on file for a single k-point. Supports plain Fortran IO and netcdf with ETSF-IO format
INPUTS
kss_unt=The unit number of the file ikpt=The index of the k-point kpoint(3)=The k-point in reduced coordinates. kss_npw=Number of planewaves used for the wavefunctions in the KSS files. npw_k=Number of planewaves at this k-point in the k-centered basis set used in abinit (ecut). ecut=cutoff energy used in abinit. rprimd(3,3)=dimensional primitive translations for real space (bohr). Psps<Pseudopotential_type>=Datatype gathering data on the Pseudopotentials. iomode=Input variables specifying the fileformat. (0-->Fortran,3-->netcdf with ETSF-IO). gbig(3,kss_npw)=Set of G-vectors used in the KSS file.
OUTPUT
Only writing.
SOURCE
352 subroutine write_vkb(kss_unt,ikpt,kpoint,kss_npw,gbig,rprimd,Psps,iomode) 353 354 !Arguments ------------------------------------ 355 !scalars 356 integer,intent(in) :: ikpt,iomode,kss_npw,kss_unt 357 type(Pseudopotential_type),intent(in) :: Psps 358 !arrays 359 integer,intent(in) :: gbig(3,kss_npw) 360 real(dp),intent(in) :: kpoint(3),rprimd(3,3) 361 362 !Local variables------------------------------- 363 !scalars 364 integer :: itypat,il,ig,mpsang,ntypat 365 !array 366 real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:) 367 real(dp),allocatable :: dum_vkbsign(:,:) 368 #ifdef HAVE_NETCDF 369 integer :: ncerr,varid 370 real(dp),allocatable,target :: vkb_tgt(:,:,:,:), vkbd_tgt(:,:,:,:) 371 #endif 372 373 ! ********************************************************************* 374 375 mpsang = Psps%mpsang; ntypat = Psps%ntypat 376 377 ABI_MALLOC(vkb ,(kss_npw,ntypat,mpsang)) 378 ABI_MALLOC(vkbd,(kss_npw,ntypat,mpsang)) 379 ABI_MALLOC(dum_vkbsign,(ntypat,mpsang)) 380 381 call kss_calc_vkb(Psps,kpoint,kss_npw,gbig,rprimd,dum_vkbsign,vkb,vkbd) 382 ABI_FREE(dum_vkbsign) 383 384 SELECT CASE (iomode) 385 386 CASE (IO_MODE_FORTRAN) 387 do itypat=1,ntypat 388 do il=1,mpsang 389 write(kss_unt) (vkb (ig,itypat,il),ig=1,kss_npw) 390 write(kss_unt) (vkbd(ig,itypat,il),ig=1,kss_npw) 391 end do 392 end do 393 394 #ifdef HAVE_NETCDF 395 CASE (IO_MODE_ETSF) 396 ABI_MALLOC(vkb_tgt ,(kss_npw,1,mpsang,ntypat)) 397 ABI_MALLOC(vkbd_tgt,(kss_npw,1,mpsang,ntypat)) 398 do itypat=1,ntypat 399 do il=1,mpsang 400 do ig=1,kss_npw 401 vkb_tgt (ig,1,il,itypat)=vkb (ig,itypat,il) 402 vkbd_tgt(ig,1,il,itypat)=vkbd(ig,itypat,il) 403 end do 404 end do 405 end do 406 407 ! FIXME: Multiple projectors 408 !ABI_MALLOC(vkbd, (npw, psps%lnmax, cryst%ntypat)) 409 !call calc_vkb(cryst,psps,kpoint,npw,npw,gvec,vkbsign,vkb,vkbd) 410 !ABI_FREE(vkbsign) 411 !ABI_FREE(vkb) 412 !ABI_FREE(vkbd) 413 414 ! The shape of the variable on disk is (Fortran API): 415 ! (max_number_of_coefficients, number_of_kpoints, max_number_of_projectors, 416 ! max_number_of_angular_momenta, number_of_atom_species) 417 varid = nctk_idname(kss_unt, "kb_formfactors") 418 ncerr = nf90_put_var(kss_unt, varid, vkb_tgt, start=[1,ikpt,1,1,1], count=[kss_npw,1,1,mpsang,ntypat]) 419 NCF_CHECK(ncerr) 420 421 varid = nctk_idname(kss_unt, "kb_formfactor_derivative") 422 ncerr = nf90_put_var(kss_unt, varid, vkbd_tgt, start=[1,ikpt,1,1,1], count=[kss_npw,1,1,mpsang,ntypat]) 423 NCF_CHECK(ncerr) 424 425 ABI_FREE(vkb_tgt) 426 ABI_FREE(vkbd_tgt) 427 #endif 428 429 CASE DEFAULT 430 ABI_ERROR(sjoin("Unsupported value for iomode:", itoa(iomode))) 431 END SELECT 432 433 ABI_FREE(vkb) 434 ABI_FREE(vkbd) 435 436 end subroutine write_vkb