TABLE OF CONTENTS
ABINIT/m_positron [ Modules ]
NAME
m_positron
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (GJ, MT, JW) 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_positron 23 24 use defs_basis 25 use m_efield 26 use m_errors 27 use m_abicore 28 use m_energies 29 use m_wffile 30 use m_electronpositron 31 use m_hdr 32 use m_xmpi 33 use m_bandfft_kpt 34 use m_dtset 35 use m_dtfil 36 use m_extfpmd 37 38 use defs_datatypes, only : pseudopotential_type 39 use defs_abitypes, only : MPI_type 40 use m_special_funcs, only : sbf8 41 use m_ioarr, only : ioarr, read_rhor 42 use m_pawang, only : pawang_type 43 use m_paw_sphharm, only : realgaunt 44 use m_pawrad, only : pawrad_type, simp_gen, nderiv_gen 45 use m_pawtab, only : pawtab_type 46 use m_paw_ij, only : paw_ij_type 47 use m_pawfgrtab,only : pawfgrtab_type 48 use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_alloc, pawrhoij_free,& 49 pawrhoij_nullify, pawrhoij_gather, pawrhoij_inquire_dim, pawrhoij_symrhoij 50 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_send, & 51 pawcprj_mpi_recv, pawcprj_free, pawcprj_copy, pawcprj_bcast 52 use m_pawfgr, only : pawfgr_type 53 use m_paw_nhat, only : pawmknhat 54 use m_fock, only : fock_type 55 use m_kg, only : getcut 56 use defs_wvltypes, only : wvl_data 57 use m_spacepar, only : hartre 58 use m_mkrho, only : initro 59 use m_paw_occupancies, only : initrhoij, pawaccrhoij 60 use m_gammapositron, only : gammapositron, gammapositron_fft 61 use m_forstr, only : forstr 62 use m_pawxc, only : pawxcsum 63 use m_paw_denpot, only : pawdensities 64 use m_drivexc, only : mkdenpos 65 66 use m_paw_sphharm, only : initylmr 67 use m_pawpsp, only : pawpsp_read_corewf 68 use m_crystal, only : crystal_t 69 use m_mpinfo, only : ptabs_fourdp,set_mpi_enreg_fft,unset_mpi_enreg_fft,destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle 70 use m_io_tools, only : open_file,close_unit,get_unit 71 use m_fftcore, only : sphereboundary 72 use m_prep_kgb, only : prep_fourwf 73 use m_fft, only : fourwf, fourdp 74 use m_cgprj, only : ctocprj 75 76 implicit none 77 78 private
ABINIT/posdoppler [ Functions ]
NAME
posdoppler
FUNCTION
Calculate the momentum distribution annihilating electrons-positron (Doppler broadening)
INPUTS
cg(2,mcg)=planewave coefficients of wavefunctions.
cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk>
and each |p_lmn> non-local projector
Cryst<Crystal_structure> = Info on unit cell and its symmetries
dimcprj(natom)=array of dimensions of array cprj (not ordered)
dtfil <type(datafiles_type)>=variables related to files
| unpaw=unit number for temporary PAW files
dtset <type(dataset_type)>=all input variables for this dataset
| istwfk=input option=1 parameter that describes the storage of wfs
| mband=maximum number of bands
| mgfft=maximum size of 1D FFTs for the "coarse" grid
| mkmem=number of k points treated by this node.
| mpw=maximum dimensioned size of npw
| natom=number of atoms
| nband=number of bands at each k point
| ngfft=contain all needed information about 3D FFT (coarse grid)
| nkpt=number of k points
| nspden=number of spin-density components
| nspinor=number of spinorial components of the wavefunctions
| nsppol=1 for unpolarized, 2 for spin-polarized
| usepaw=flag for PAW
| gpu_option=GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
| wtk(=weights associated with various k points
filpsp(ntypat)=name(s) of the pseudopotential file(s)
kg(3,mpw*mkmem)=reduced planewave coordinates.
mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
mpi_enreg= information about MPI parallelization
my_natom=number of atoms treated by current processor
n3xccc= dimension of the xccc3d array (0 or nfft).
nfft= number of FFT grid points
ngfft(18)= contain all needed information about 3D FFT
nhat(nfft,nspden)=charge compensation density (content depends on electronpositron%particle)
npwarr(nkpt)=number of planewaves in basis at this k point
occ(mband*nkpt*nsppol)=occupancies for each band and k point
pawang <type(pawang)>=paw angular mesh and related data
pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
rhor(nfft,nspden)=total electron/positron density (content depends on electronpositron%particle)
xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
OUTPUT
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
TODO
print a warning if the core wave function is not localized in the PAW sphere implement PAW on-site contribution for state-independent scheme
SOURCE
1833 !Macro to go from row-column indexing to combined indexing 1834 #define RCC(glmn,hlmn) max(glmn,hlmn)*(max(glmn,hlmn)-1)/2+min(glmn,hlmn) 1835 !Macro to go from l,m angular momentum indexing to combined indexing 1836 #define LMC(lval,mval) lval*lval+lval+mval+1 1837 1838 subroutine posdoppler(cg,cprj,Crystal,dimcprj,dtfil,dtset,electronpositron,& 1839 & filpsp,kg,mcg,mcprj,mpi_enreg,my_natom,& 1840 & n3xccc,nfft,ngfft,nhat,npwarr,occ,pawang,pawrad,& 1841 & pawrhoij,pawtab,rhor,xccc3d) 1842 1843 !Arguments ------------------------------------ 1844 !scalars 1845 integer,intent(in) :: mcg,mcprj,my_natom,n3xccc,nfft 1846 type(crystal_t) :: Crystal 1847 type(datafiles_type),intent(in) :: dtfil 1848 type(dataset_type),intent(in) :: dtset 1849 type(electronpositron_type),pointer :: electronpositron 1850 type(MPI_type),intent(inout) :: mpi_enreg 1851 type(pawang_type),intent(in) :: pawang 1852 !arrays 1853 integer,intent(in) :: dimcprj(dtset%natom) 1854 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),ngfft(18),npwarr(dtset%nkpt) 1855 real(dp),intent(in) :: nhat(nfft,dtset%nspden*dtset%usepaw),xccc3d(n3xccc) 1856 real(dp),intent(in),target :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 1857 real(dp),intent(in),target :: rhor(nfft,dtset%nspden) 1858 real(dp),intent(inout),target :: cg(2,mcg) 1859 character(len=fnlen),intent(in) :: filpsp(dtset%ntypat) 1860 type(pawcprj_type),target :: cprj(dtset%natom,mcprj) 1861 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw) 1862 type(pawrhoij_type),intent(in),target :: pawrhoij(my_natom*dtset%usepaw) 1863 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 1864 1865 !Local variables------------------------------- 1866 !scalars 1867 integer :: accessfil,basis_size,bandpp,bdtot_index,bdtot_index_pos,blocksize,cplex,cplex_rhoij 1868 integer :: glmij,i0lmn,i1,i2,i3,iat,iatm,iatom 1869 integer :: ib,ib_cprj,ib_cprj_pos,ib_pos,ibg,ibg_pos 1870 integer :: iblock,iblock_pos,ibpp,ibpp_pos 1871 integer :: icg,icg_pos,id1,id2,id3,ierr,ig1,ig2,ig3,igamma,ii,ikg,ikg_pos,ikpt 1872 integer :: ikpt_pos,il,ilm,ilmn,iln,indx,indx0,iorder_cprj,iproc,ir,isppol,isppol_pos,istwf_k 1873 integer :: istwf_k_pos,itypat,iwarn,iwavef,iwavef_pos,j2,j3,jj,jkpt,jl,jlm,jlmn,jln 1874 integer :: klm,kln,klmn,l_size,l_size_max,ll,llmax,llmin,lm,lmn_size,lmn_size_c,lmn2_size 1875 integer :: mband_cprj,mband_cprj_pos,mcg_pos 1876 integer :: mcprj_k,mcprj_k_pos,me_band,me_fft,me_kpt,me_kptband 1877 integer :: mesh_size,meshsz,mm,my_ngrid,my_nspinor,my_nsppol,my_n2,n1,n2,n3,n4,n5,n6 1878 integer :: nband_cprj_eff_pos,nband_cprj_k,nband_cprj_k_pos 1879 integer :: nband_eff_pos,nband_k,nband_k_pos 1880 integer :: nblock_band,nblock_band_eff_pos,nkpt 1881 integer :: nproc_band,nproc_fft,nproc_spkpt,nproc_kptband,npw_k,npw_k_pos 1882 integer :: nspden_rhoij,option,tag,unit_doppler 1883 integer :: tim_fourdp=0,tim_fourwf=-36 1884 integer :: ylmr_normchoice,ylmr_npts,ylmr_option 1885 logical,parameter :: include_nhat_in_gamma=.false.,state_dependent=.true. 1886 logical,parameter :: kgamma_only_positron=.true.,wf_conjugate=.false. 1887 logical :: cprj_paral_band,ex,mykpt,mykpt_pos,use_timerev,use_zeromag,abinitcorewf,xmlcorewf 1888 real(dp) :: arg,bessarg,cpi,cpr,cp11,cp12,cp21,cp22,gammastate,intg 1889 real(dp) :: lambda_v1,lambda_v2,lambda_core,lambda_pw,occ_el,occ_pos 1890 real(dp) :: pnorm,pr,rate,rate_ipm,ratec,ratec_ipm,rate_paw,rate_paw_ipm 1891 real(dp) :: scale_,units_,weight,weight_pos,wf_fact,wtk_k,wtk_k_pos,vec 1892 character(len=fnlen) :: filename,filename_dop 1893 character(len=500) :: msg 1894 type(bandfft_kpt_type),pointer :: bandfft_kpt_el,bandfft_kpt_pos 1895 type(MPI_type) :: mpi_enreg_seq 1896 type(wffile_type) :: wff 1897 !arrays 1898 integer,allocatable :: gbound(:,:),gbound_pos(:,:),kg_k(:,:),kg_k_pos(:,:) 1899 integer,allocatable :: lcor(:),lmncmax(:),my_ffttab(:),my_gridtab(:),ncor(:),nphicor(:) 1900 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1901 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1902 logical,allocatable :: have_intc(:,:,:),have_rad(:,:) 1903 real(dp) :: buf(4),contrib(2),cp(2),cp_pos(2),expipr(2),pbn(3),pcart(3) 1904 real(dp) :: radsumnfftc(2),ylmgr(1,1,0),ylmr_nrm(1) 1905 real(dp),allocatable :: cwaveg(:,:),cwaveg_pos(:,:),cwaver(:),cwaver_pos(:),cwaver_pos_block(:) 1906 real(dp),allocatable :: cg_k_pos(:,:),cwaveaug(:,:,:,:),cwaveaug_pos(:,:,:,:) 1907 real(dp),allocatable :: denpot_dum(:,:,:),energycor(:),ff(:),fofgout_dum(:,:) 1908 real(dp),allocatable :: gamma(:,:),intc(:,:,:),j_bessel(:,:),jbes(:),mpibuf(:,:) 1909 real(dp),allocatable :: occ_k(:),occ_k_pos(:),pcart_k(:,:) 1910 real(dp),allocatable :: radint1(:,:),radint2(:,:),radint3(:,:) 1911 real(dp),allocatable :: radsumnfft1(:,:),radsumnfft2(:,:),radsumnfft3(:,:) 1912 real(dp),allocatable :: rho_contrib(:),rho_contrib_g(:,:) 1913 real(dp),allocatable :: rho_contrib_paw1(:,:),rho_contrib_paw2(:,:),rho_contrib_paw3(:,:) 1914 real(dp),allocatable :: rho_moment_v1(:,:),rho_moment_v2(:,:) 1915 real(dp),allocatable :: rho_moment_core(:,:),rho_moment_k(:),rho_moment_k2(:) 1916 real(dp),allocatable :: rho_pw(:,:),rhor_dop_el(:) 1917 real(dp),allocatable :: rhocorej(:),rhoe(:,:),rhop(:,:),ylmp(:) 1918 real(dp),pointer :: cg_pos_ptr(:,:),cg_ptr(:,:),occ_ptr(:),occ_pos_ptr(:) 1919 real(dp),pointer :: rhor_(:,:),rhor_ep_(:,:) 1920 complex(dpc) :: ifac ! (-i)^L mod 4 1921 complex(dpc),dimension(0:3) :: ilfac(0:3)=(/(1.0,0.0),(0.0,-1.0),(-1.0,0.0),(0.0,1.0)/) 1922 type(coeff1_type),allocatable :: gammastate_c(:) 1923 type(coeffi2_type),allocatable :: indlmncor(:) 1924 type(coeff2_type),allocatable :: phicor(:) 1925 type(coeff6_type),allocatable :: radsum1(:),radsum2(:),radsum3(:) 1926 type(coeff7_type),allocatable :: radsumc(:) 1927 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_k_pos(:,:),cprj_pos(:,:) 1928 type(pawcprj_type),pointer :: cprj_pos_ptr(:,:),cprj_ptr(:,:) 1929 type(pawrhoij_type),allocatable :: pawrhoij_dop_el(:) 1930 type(pawrhoij_type),pointer :: pawrhoij_ptr(:),pawrhoij_all(:),pawrhoij_ep_all(:) 1931 1932 ! ************************************************************************* 1933 1934 DBG_ENTER("COLL") 1935 1936 !Compatibility tests 1937 if (.not.associated(electronpositron)) then 1938 msg='electronpositron variable must be associated!' 1939 ABI_BUG(msg) 1940 end if 1941 if (allocated(mpi_enreg%proc_distrb)) then 1942 do isppol=1,dtset%nsppol 1943 do ikpt=1,dtset%nkpt 1944 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 1945 if (any(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)/=mpi_enreg%proc_distrb(ikpt,1,isppol))) then 1946 msg='proc_distrib cannot be distributed over bands!' 1947 ABI_BUG(msg) 1948 end if 1949 end do 1950 end do 1951 end if 1952 if (dtset%nspinor==2) then 1953 msg='Doppler broadening not available for spinorial wave functions (nspinor=2)!' 1954 ABI_BUG(msg) 1955 end if 1956 if (mcprj==0) then 1957 msg='<p|Psi> (cprj) datastructure must be kept in memory (see pawusecp input keyword)!' 1958 ABI_BUG(msg) 1959 end if 1960 if (dtset%usepaw==0) then 1961 write(msg,'(5a)') 'Momentum distribution of annihilating electron-positron pairs',ch10,& 1962 & 'in the Norm-conserving Pseudopotential formalism is incomplete!',ch10,& 1963 & 'No core contribution is included.' 1964 ABI_WARNING(msg) 1965 end if 1966 if (any(dtset%nband(:)/=dtset%nband(1))) then 1967 write(msg,'(a)') 'Number of bands has to be the same for all k-points!' 1968 ABI_BUG(msg) 1969 end if 1970 if (dtset%usepaw==1) then 1971 if (size(pawrhoij)/=mpi_enreg%my_natom) then 1972 write(msg,'(a)') 'wrong size for pawrhoij! ' 1973 ABI_BUG(msg) 1974 end if 1975 end if 1976 1977 !Various initializations 1978 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1979 n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6) 1980 id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2 1981 iorder_cprj=0 ; cplex=2 ; iwarn=1 1982 wf_fact=one;if (wf_conjugate) wf_fact=-one 1983 nkpt=dtset%nkpt 1984 1985 !Manage kpt/spin parallelism 1986 ABI_MALLOC(my_gridtab,(nkpt)) 1987 my_gridtab=0 1988 do ii=1,nkpt 1989 if (any(mpi_enreg%my_isppoltab(:)==1)) my_gridtab(ii)=mpi_enreg%my_kpttab(ii) 1990 end do 1991 my_ngrid=count(my_gridtab(:)/=0) 1992 my_nsppol=sum(mpi_enreg%my_isppoltab(:)) 1993 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 1994 1995 !Parallel settings 1996 if (mpi_enreg%paral_kgb/=0) then 1997 nproc_spkpt=mpi_enreg%nproc_spkpt 1998 nproc_band=mpi_enreg%nproc_band 1999 nproc_fft=mpi_enreg%nproc_fft 2000 nproc_kptband=xmpi_comm_size(mpi_enreg%comm_kptband) 2001 me_kpt=mpi_enreg%me_kpt 2002 me_band=mpi_enreg%me_band 2003 me_fft=mpi_enreg%me_fft 2004 me_kptband=xmpi_comm_rank(mpi_enreg%comm_kptband) 2005 bandpp=mpi_enreg%bandpp 2006 my_n2=n2/nproc_fft 2007 accessfil=IO_MODE_FORTRAN;if(nproc_fft>1)accessfil=IO_MODE_MPI 2008 else 2009 nproc_spkpt=mpi_enreg%nproc_spkpt 2010 nproc_band=1;nproc_fft=1 2011 nproc_kptband=nproc_spkpt 2012 me_band=0;me_fft=0 2013 me_kpt=mpi_enreg%me_kpt 2014 me_kptband=me_kpt 2015 bandpp=1 ; my_n2=n2 2016 accessfil=IO_MODE_FORTRAN 2017 end if 2018 blocksize=nproc_band*bandpp 2019 nblock_band=dtset%nband(1)/blocksize 2020 2021 !Select density according to nhat choice0 2022 if (dtset%usepaw==0.or.include_nhat_in_gamma) then 2023 rhor_ => rhor 2024 rhor_ep_ => electronpositron%rhor_ep 2025 else 2026 ABI_MALLOC(rhor_,(nfft,dtset%nspden)) 2027 ABI_MALLOC(rhor_ep_,(nfft,dtset%nspden)) 2028 rhor_=rhor-nhat 2029 rhor_ep_=electronpositron%rhor_ep-electronpositron%nhat_ep 2030 end if 2031 2032 !Select type(s) of enhancement factor 2033 igamma=0 2034 if (electronpositron%ixcpositron==-1) igamma=0 2035 if (electronpositron%ixcpositron== 1) igamma=2 2036 if (electronpositron%ixcpositron== 2) igamma=4 2037 if (electronpositron%ixcpositron== 3) igamma=2 2038 if (electronpositron%ixcpositron==11) igamma=3 2039 if (electronpositron%ixcpositron==31) igamma=3 2040 2041 !Select electronic and positronic states 2042 if (electronpositron%particle==EP_ELECTRON) then !we should not be in this case 2043 cg_ptr => electronpositron%cg_ep 2044 cprj_ptr => electronpositron%cprj_ep 2045 occ_ptr => electronpositron%occ_ep 2046 pawrhoij_ptr => electronpositron%pawrhoij_ep 2047 cg_pos_ptr => cg 2048 cprj_pos_ptr => cprj 2049 occ_pos_ptr => occ 2050 end if 2051 if (electronpositron%particle==EP_POSITRON) then 2052 cg_ptr => cg 2053 cprj_ptr => cprj 2054 occ_ptr => occ 2055 pawrhoij_ptr => pawrhoij 2056 cg_pos_ptr => electronpositron%cg_ep 2057 cprj_pos_ptr => electronpositron%cprj_ep 2058 occ_pos_ptr => electronpositron%occ_ep 2059 end if 2060 2061 !Determine if cprj datastructures are distributed over bands 2062 mband_cprj=size(cprj_ptr,2)/(my_nspinor*dtset%mkmem*dtset%nsppol) 2063 mband_cprj_pos=size(cprj_pos_ptr,2)/(my_nspinor*dtset%mkmem*dtset%nsppol) 2064 cprj_paral_band=(mband_cprj<dtset%mband) 2065 2066 !Get the distrib associated with the fft_grid 2067 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 2068 2069 !=============================================================================== 2070 !================ Calculate the PAW on-site constants ========================== 2071 2072 if (dtset%usepaw==1) then 2073 2074 ylmr_normchoice = 0 ! input to initylmr are normalized 2075 ylmr_npts = 1 ! only 1 point to compute in initylmr 2076 ylmr_nrm(1) = one ! weight of normed point for initylmr 2077 ylmr_option = 1 ! compute only ylm's in initylmr 2078 2079 !Prepare radial integral for PAW correction for each atom type 2080 ABI_MALLOC(radsum1,(dtset%ntypat)) 2081 ABI_MALLOC(radsum2,(dtset%ntypat)) 2082 ABI_MALLOC(radsum3,(dtset%ntypat)) 2083 ABI_MALLOC(radsumc,(dtset%ntypat)) 2084 2085 ABI_MALLOC(indlmncor,(dtset%ntypat)) 2086 ABI_MALLOC(phicor,(dtset%ntypat)) 2087 ABI_MALLOC(gammastate_c,(dtset%natom)) 2088 ABI_MALLOC(nphicor,(dtset%ntypat)) 2089 ABI_MALLOC(lmncmax,(dtset%ntypat)) 2090 2091 ! Reading of core wave functions 2092 if (mpi_enreg%me_cell==0) then 2093 do itypat=1,dtset%ntypat 2094 filename=trim(filpsp(itypat)) ; iln=len(trim(filename)) 2095 abinitcorewf=.false. ; if (iln>3) abinitcorewf=(filename(iln-6:iln)=='.abinit') 2096 xmlcorewf=.false. ; if (iln>3) xmlcorewf=(filename(iln-3:iln)=='.xml') 2097 if ((.not.xmlcorewf).and.(.not.abinitcorewf)) filename=filename(1:iln)//'.corewf' 2098 if (abinitcorewf) filename=filename(1:iln-6)//'corewf.abinit' 2099 if (xmlcorewf) filename=filename(1:iln-3)//'corewf.xml' 2100 inquire(file=filename,exist=ex) 2101 if (.not.ex) then 2102 write(unit=filename,fmt='(a,i1)') 'corewf.abinit',itypat 2103 inquire(file=filename,exist=ex) 2104 end if 2105 if (.not.ex) then 2106 write(msg,'(3a)') 'Core wave-functions file is missing!',ch10,& 2107 & 'Looking for: psp-name.corewf[.xml][.abinit] or corewf.dat' 2108 ABI_ERROR(msg) 2109 end if 2110 call pawpsp_read_corewf(energycor,indlmncor(itypat)%value,lcor,lmncmax(itypat),& 2111 & ncor,nphicor(itypat),pawrad(itypat),phicor(itypat)%value,& 2112 & filename=filename) 2113 ! The following arrays are not used anymore 2114 if (nphicor(itypat)>0) then 2115 ABI_FREE(energycor) 2116 ABI_FREE(lcor) 2117 ABI_FREE(ncor) 2118 end if 2119 end do 2120 end if 2121 if (mpi_enreg%nproc_cell>1) then 2122 call xmpi_bcast(indlmncor,0,mpi_enreg%comm_cell,ierr) 2123 call xmpi_bcast(phicor,0,mpi_enreg%comm_cell,ierr) 2124 call xmpi_bcast(nphicor,0,mpi_enreg%comm_cell,ierr) 2125 call xmpi_bcast(lmncmax,0,mpi_enreg%comm_cell,ierr) 2126 end if 2127 2128 do itypat=1,dtset%ntypat 2129 2130 mesh_size = pawtab(itypat)%mesh_size 2131 l_size = pawtab(itypat)%l_size 2132 lmn_size = pawtab(itypat)%lmn_size 2133 lmn2_size = pawtab(itypat)%lmn2_size 2134 basis_size = pawtab(itypat)%basis_size 2135 2136 lmn_size_c=lmncmax(itypat) 2137 llmax=maxval(indlmncor(itypat)%value(1,1:lmn_size_c)) 2138 l_size_max=max(l_size,2*llmax+1) 2139 2140 ABI_MALLOC(j_bessel,(mesh_size,l_size_max)) 2141 ABI_MALLOC(ylmp,(l_size_max*l_size_max)) 2142 ABI_MALLOC(have_intc,(l_size_max,basis_size,nphicor(itypat))) 2143 ABI_MALLOC(intc,(l_size_max,basis_size,nphicor(itypat))) 2144 ABI_MALLOC(have_rad,(l_size,pawtab(itypat)%ij_size)) 2145 ABI_MALLOC(radint1,(l_size,pawtab(itypat)%ij_size)) 2146 ABI_MALLOC(radint2,(l_size,pawtab(itypat)%ij_size)) 2147 ABI_MALLOC(radint3,(l_size,pawtab(itypat)%ij_size)) 2148 2149 ABI_MALLOC(radsumc(itypat)%value,(2,lmn_size,lmn_size_c,n1,my_n2,n3,my_ngrid)) 2150 ABI_MALLOC(radsum1(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid)) 2151 ABI_MALLOC(radsum2(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid)) 2152 ABI_MALLOC(radsum3(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid)) 2153 radsumc(itypat)%value=zero 2154 radsum1(itypat)%value=zero 2155 radsum2(itypat)%value=zero 2156 radsum3(itypat)%value=zero 2157 2158 ABI_MALLOC(jbes,(l_size_max)) 2159 ABI_MALLOC(ff,(mesh_size)) 2160 meshsz=pawrad(itypat)%int_meshsz 2161 if (meshsz>mesh_size) ff(meshsz+1:mesh_size)=zero 2162 2163 indx=0;jkpt=0 2164 do ikpt=1,nkpt 2165 if (my_gridtab(ikpt)==0) cycle 2166 jkpt=jkpt+1 2167 do i3=1,n3 2168 ig3=i3-(i3/id3)*n3-1 2169 do i2=1,n2 2170 if (me_fft/=fftn2_distrib(i2)) cycle 2171 j2=ffti2_local(i2) 2172 indx=n1*(my_n2*(i3-1)+(j2-1)) 2173 ig2=i2-(i2/id2)*n2-1 2174 do i1=1,n1 2175 ig1=i1-(i1/id1)*n1-1 2176 indx=indx+1;if (mod(indx-1,nproc_band)/=me_band) cycle 2177 2178 pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+& 2179 & Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+& 2180 & Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt)) 2181 pnorm=dsqrt(dot_product(pcart,pcart)) 2182 2183 if (pnorm < tol12) then 2184 pbn(:) = zero 2185 ylmp(:) = zero 2186 ylmp(1) = 1.d0/sqrt(four_pi) 2187 else 2188 pbn(:) = pcart(:)/pnorm ! unit vector 2189 call initylmr(l_size,ylmr_normchoice,ylmr_npts,ylmr_nrm,ylmr_option,pbn(:),ylmp(:),ylmgr) 2190 end if 2191 2192 pnorm=two_pi*pnorm ! re-normed for call to bessel 2193 do ir = 1, mesh_size 2194 bessarg = pnorm*pawrad(itypat)%rad(ir) 2195 call sbf8(l_size,bessarg,jbes) 2196 j_bessel(ir,:)=jbes(:) 2197 end do 2198 2199 ! ===== Core part ===== 2200 ! Need intc=\int phi phi_core jl (pr) dr 2201 2202 have_intc(:,:,:)=.FALSE. ; intc(:,:,:)=zero 2203 2204 do jlmn = 1,lmn_size_c 2205 jln = indlmncor(itypat)%value(5,jlmn) 2206 jlm = indlmncor(itypat)%value(4,jlmn) 2207 jl = indlmncor(itypat)%value(1,jlmn) 2208 do ilmn = 1,lmn_size 2209 iln = pawtab(itypat)%indlmn(5,ilmn) 2210 ilm = pawtab(itypat)%indlmn(4,ilmn) 2211 il = pawtab(itypat)%indlmn(1,ilmn) 2212 2213 llmin = abs(il-jl) 2214 llmax = il+jl 2215 klm = RCC(ilm,jlm) 2216 do ll=llmin,llmax,2 2217 ifac=ilfac(mod(ll,4)) 2218 2219 if (.not.have_intc(ll+1,iln,jln)) then 2220 ff(1:mesh_size)=(pawtab(itypat)%phi(1:mesh_size,iln)*phicor(itypat)%value(1:mesh_size,jln))& 2221 & *j_bessel(1:mesh_size,ll+1) 2222 call simp_gen(intg,ff,pawrad(itypat)) 2223 intc(ll+1,iln,jln)=intg 2224 have_intc(ll+1,iln,jln)=.true. 2225 end if 2226 2227 do mm=-ll,ll 2228 lm = LMC(ll,mm) 2229 glmij=pawang%gntselect(lm,klm) 2230 if (glmij>0) then 2231 arg=ylmp(lm)*pawang%realgnt(glmij)*intc(ll+1,iln,jln) 2232 radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt) = & 2233 & radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt)+arg*real(ifac) 2234 radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt) = & 2235 & radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt)+arg*aimag(ifac) 2236 end if 2237 end do !mm 2238 end do !ll 2239 end do !ilmn 2240 end do !jlmn 2241 2242 ! ===== Valence part ===== 2243 ! Need int1=\int phi_i phi_j jl (pr) dr 2244 ! and int2=\int tphi_i tphi_j jl (pr) dr 2245 2246 have_rad(:,:)= .FALSE.;radint1=zero;radint2=zero;radint3=zero 2247 2248 do klmn=1,pawtab(itypat)%lmn2_size 2249 klm=pawtab(itypat)%indklmn(1,klmn);kln=pawtab(itypat)%indklmn(2,klmn) 2250 llmin=pawtab(itypat)%indklmn(3,klmn);llmax=pawtab(itypat)%indklmn(4,klmn) 2251 2252 do ll=llmin,llmax,2 2253 ifac=ilfac(mod(ll,4)) 2254 2255 if (.not.have_rad(ll+1,kln)) then 2256 ff(1:mesh_size)=pawtab(itypat)%phiphj(1:mesh_size,kln)*j_bessel(1:mesh_size,ll+1) 2257 call simp_gen(intg,ff,pawrad(itypat)) 2258 radint1(ll+1,kln)=intg 2259 ff(1:mesh_size)=pawtab(itypat)%tphitphj(1:mesh_size,kln)*j_bessel(1:mesh_size,ll+1) 2260 call simp_gen(intg,ff,pawrad(itypat)) 2261 radint2(ll+1,kln)=intg 2262 ff(1:mesh_size)=(pawtab(itypat)%phiphj (1:mesh_size,kln) & 2263 & -pawtab(itypat)%tphitphj(1:mesh_size,kln))& 2264 & *j_bessel(1:mesh_size,ll+1) 2265 call simp_gen(intg,ff,pawrad(itypat)) 2266 radint3(ll+1,kln)=intg 2267 have_rad(ll+1,kln)=.true. 2268 end if 2269 2270 do mm=-ll,ll 2271 lm = LMC(ll,mm) 2272 glmij=pawang%gntselect(lm,klm) 2273 if (glmij>0) then 2274 arg=ylmp(lm)*pawang%realgnt(glmij) 2275 radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt) = & 2276 & radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint1(ll+1,kln) 2277 radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt) = & 2278 & radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint1(ll+1,kln) 2279 radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt) = & 2280 & radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint2(ll+1,kln) 2281 radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt) = & 2282 & radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint2(ll+1,kln) 2283 radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt) = & 2284 & radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint3(ll+1,kln) 2285 radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt) = & 2286 & radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint3(ll+1,kln) 2287 end if 2288 end do !mm 2289 end do !ll 2290 end do !klmn 2291 2292 end do ! end loop over i1 2293 end do ! end loop over i2 2294 end do ! end loop over i3 2295 end do ! end loop over ikpt 2296 2297 ABI_FREE(ff) 2298 ABI_FREE(jbes) 2299 2300 ABI_FREE(j_bessel) 2301 ABI_FREE(ylmp) 2302 2303 ABI_FREE(intc) 2304 ABI_FREE(have_intc) 2305 2306 ABI_FREE(radint1) 2307 ABI_FREE(radint2) 2308 ABI_FREE(radint3) 2309 ABI_FREE(have_rad) 2310 2311 call xmpi_sum(radsumc(itypat)%value,mpi_enreg%comm_band,ierr) 2312 call xmpi_sum(radsum1(itypat)%value,mpi_enreg%comm_band,ierr) 2313 call xmpi_sum(radsum2(itypat)%value,mpi_enreg%comm_band,ierr) 2314 call xmpi_sum(radsum3(itypat)%value,mpi_enreg%comm_band,ierr) 2315 2316 end do ! end loop over atom types 2317 end if ! PAW 2318 2319 !Allocate main memory 2320 ABI_MALLOC(rho_contrib,(cplex*nfft)) 2321 ABI_MALLOC(rho_contrib_g,(cplex,nfft)) 2322 ABI_MALLOC(rho_contrib_paw1,(cplex,nfft)) 2323 ABI_MALLOC(rho_contrib_paw2,(cplex,nfft)) 2324 ABI_MALLOC(rho_contrib_paw3,(cplex,nfft)) 2325 2326 ABI_MALLOC(rho_moment_v1,(nfft,my_ngrid)) 2327 ABI_MALLOC(rho_moment_v2,(nfft,my_ngrid)) 2328 ABI_MALLOC(rho_moment_core,(nfft,my_ngrid)) 2329 ABI_MALLOC(rho_pw,(nfft,my_ngrid)) 2330 rho_moment_v1=zero;rho_moment_v2=zero 2331 rho_pw=zero;rho_moment_core=zero 2332 2333 !Prepare gamma(r) for the state independent scheme 2334 ABI_MALLOC(gamma,(nfft,2)) 2335 if (.not.state_dependent) then 2336 ABI_MALLOC(rhoe,(nfft,1)) 2337 ABI_MALLOC(rhop,(nfft,1)) 2338 if (electronpositron%particle==EP_ELECTRON) then 2339 rhoe(:,1)=rhor_ep_(:,1);rhop(:,1)=rhor_(:,1) 2340 else if (electronpositron%particle==EP_POSITRON) then 2341 rhoe(:,1)=rhor_(:,1);rhop(:,1)=rhor_ep_(:,1) 2342 end if 2343 call mkdenpos(iwarn,nfft,1,1,rhoe(:,1),dtset%xc_denpos) 2344 call mkdenpos(iwarn,nfft,1,1,rhop(:,1),dtset%xc_denpos) 2345 call gammapositron_fft(electronpositron,gamma,Crystal%gprimd,igamma,mpi_enreg,& 2346 & n3xccc,nfft,ngfft,rhoe(:,1),rhop(:,1),xccc3d) 2347 ABI_FREE(rhoe) 2348 ABI_FREE(rhop) 2349 else 2350 gamma=one 2351 end if 2352 2353 !Some allocations for state-dependent scheme 2354 if (state_dependent) then 2355 ! Fake MPI data to be used in poslifetime; allow only FFT parallelism 2356 call initmpi_seq(mpi_enreg_seq) 2357 mpi_enreg_seq%my_natom=dtset%natom 2358 call set_mpi_enreg_fft(mpi_enreg_seq,mpi_enreg%comm_fft,mpi_enreg%distribfft,& 2359 & mpi_enreg%me_g0,mpi_enreg%paral_kgb) 2360 ! Allocate memory for state-dependent scheme 2361 ABI_MALLOC(rhor_dop_el,(nfft)) 2362 if (dtset%usepaw==1) then 2363 ABI_MALLOC(pawrhoij_dop_el,(dtset%natom)) 2364 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,& 2365 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc) 2366 call pawrhoij_alloc(pawrhoij_dop_el,cplex_rhoij,nspden_rhoij,& 2367 dtset%nspinor,dtset%nsppol,dtset%typat,& 2368 pawtab=pawtab,use_rhoij_=1,use_rhoijp=1) 2369 ! Cancel distribution of PAW data over atomic sites 2370 ! We use here pawrhoij because polifetime routine 2371 ! detects by itself the particle described by pawrhoij 2372 if (mpi_enreg%my_natom<dtset%natom) then 2373 ABI_MALLOC(pawrhoij_all,(dtset%natom)) 2374 call pawrhoij_nullify(pawrhoij_all) 2375 call pawrhoij_gather(pawrhoij,pawrhoij_all,-1,mpi_enreg%comm_atom, & 2376 & with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.) 2377 ABI_MALLOC(pawrhoij_ep_all,(dtset%natom)) 2378 call pawrhoij_nullify(pawrhoij_ep_all) 2379 call pawrhoij_gather(electronpositron%pawrhoij_ep,pawrhoij_ep_all,-1,mpi_enreg%comm_atom, & 2380 & with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.) 2381 else 2382 pawrhoij_all => pawrhoij 2383 pawrhoij_ep_all => electronpositron%pawrhoij_ep 2384 end if 2385 end if 2386 end if 2387 2388 !============================================================================== 2389 !================ Loop over positronic states ================================= 2390 2391 !LOOP OVER k POINTS 2392 ibg_pos=0;icg_pos=0;ikg_pos=0;bdtot_index_pos=0;isppol_pos=1 2393 do ikpt_pos=1,merge(1,nkpt,kgamma_only_positron) 2394 2395 ! Extract data for this kpt_pos 2396 npw_k_pos=npwarr(ikpt_pos) 2397 wtk_k_pos=dtset%wtk(ikpt_pos); if (kgamma_only_positron) wtk_k_pos=one 2398 istwf_k_pos=dtset%istwfk(ikpt_pos) 2399 nband_k_pos=dtset%nband(ikpt_pos+(isppol_pos-1)*nkpt) 2400 nband_cprj_k_pos=nband_k_pos/nproc_band 2401 mykpt_pos=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_pos,1,nband_k_pos,& 2402 & isppol_pos,mpi_enreg%me_kpt)) 2403 2404 ! Retrieve additional data for this kpt_pos 2405 ABI_MALLOC(occ_k_pos,(nband_k_pos)) 2406 occ_k_pos(:)=occ_pos_ptr(1+bdtot_index_pos:nband_k_pos+bdtot_index_pos) 2407 nband_eff_pos=1 2408 do ib_pos=1,nband_k_pos 2409 if (occ_k_pos(ib_pos)>tol8) nband_eff_pos=ib_pos 2410 end do 2411 if (mod(nband_eff_pos,blocksize)/=0) nband_eff_pos=((nband_eff_pos/blocksize)+1)*blocksize 2412 2413 nblock_band_eff_pos=nband_eff_pos/blocksize 2414 2415 mcg_pos=npw_k_pos*my_nspinor*nband_eff_pos 2416 ABI_MALLOC(cg_k_pos,(2,mcg_pos)) 2417 2418 mcprj_k_pos=0 2419 if (dtset%usepaw==1) then 2420 nband_cprj_eff_pos=nband_eff_pos/nproc_band 2421 mcprj_k_pos=my_nspinor*nband_cprj_eff_pos 2422 ABI_MALLOC(cprj_k_pos,(dtset%natom,mcprj_k_pos)) 2423 call pawcprj_alloc(cprj_k_pos,0,dimcprj) 2424 end if 2425 2426 if (mpi_enreg%paral_kgb==0) then 2427 ABI_MALLOC(gbound_pos,(2*dtset%mgfft+8,2)) 2428 ABI_MALLOC(kg_k_pos,(3,npw_k_pos)) 2429 else if (mykpt_pos) then 2430 nullify(bandfft_kpt_pos) 2431 else 2432 ABI_MALLOC(bandfft_kpt_pos,) 2433 call bandfft_kpt_reset(bandfft_kpt_pos) 2434 end if 2435 2436 ! Exchange data (WF components) between procs 2437 if (mykpt_pos) then 2438 cg_k_pos(:,1:mcg_pos)=cg_pos_ptr(:,icg_pos+1:icg_pos+mcg_pos) 2439 if (mpi_enreg%paral_kgb==0) kg_k_pos(:,1:npw_k_pos)=kg(:,1+ikg_pos:npw_k_pos+ikg_pos) 2440 if (dtset%usepaw==1) then 2441 call pawcprj_get(Crystal%atindx1,cprj_k_pos,cprj_pos_ptr,dtset%natom,1,ibg_pos,ikpt_pos,iorder_cprj,& 2442 & isppol_pos,mband_cprj_pos,dtset%mkmem,dtset%natom,nband_cprj_eff_pos,nband_k_pos,my_nspinor,& 2443 & dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 2444 end if 2445 if (mpi_enreg%paral_kgb/=0) then 2446 jj=mpi_enreg%my_kpttab(ikpt_pos) 2447 bandfft_kpt_pos => bandfft_kpt(jj) 2448 end if 2449 do ii=0,mpi_enreg%nproc_spkpt-1 2450 if (ii/=mpi_enreg%me_kpt) then 2451 tag=ikpt_pos+(isppol_pos-1)*nkpt+2*nkpt*ii 2452 call xmpi_send(cg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr) 2453 tag=tag+nkpt*(1+2*mpi_enreg%nproc_spkpt) 2454 if (mpi_enreg%paral_kgb==0) then 2455 call xmpi_send(kg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr) 2456 else 2457 call bandfft_kpt_mpi_send(bandfft_kpt_pos,ii,tag,mpi_enreg%comm_kpt,ierr,profile='fourwf') 2458 end if 2459 if (dtset%usepaw==1) then 2460 call pawcprj_mpi_send(dtset%natom,mcprj_k_pos,dimcprj,0,cprj_k_pos,ii,mpi_enreg%comm_kpt,ierr) 2461 end if 2462 end if 2463 end do 2464 else 2465 ii=0;if (allocated(mpi_enreg%proc_distrb)) ii=mpi_enreg%proc_distrb(ikpt_pos,1,isppol_pos) 2466 tag=ikpt_pos+(isppol_pos-1)*nkpt+2*nkpt*mpi_enreg%me_kpt 2467 call xmpi_recv(cg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr) 2468 tag=tag+nkpt*(1+2*mpi_enreg%nproc_spkpt) 2469 if (mpi_enreg%paral_kgb==0) then 2470 call xmpi_recv(kg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr) 2471 else 2472 call bandfft_kpt_mpi_recv(bandfft_kpt_pos,ii,tag,mpi_enreg%comm_kpt,ierr) 2473 end if 2474 if (dtset%usepaw==1) then 2475 call pawcprj_mpi_recv(dtset%natom,mcprj_k_pos,dimcprj,0,cprj_k_pos,ii,mpi_enreg%comm_kpt,ierr) 2476 end if 2477 end if 2478 2479 if (mpi_enreg%paral_kgb==0) then 2480 call sphereboundary(gbound_pos,istwf_k_pos,kg_k_pos,dtset%mgfft,npw_k_pos) 2481 end if 2482 2483 ABI_MALLOC(cwaver_pos,(cplex*nfft)) 2484 ABI_MALLOC(cwaver_pos_block,(cplex*nfft*bandpp)) 2485 if (dtset%usepaw==1) then 2486 ABI_MALLOC(cprj_pos,(dtset%natom,my_nspinor)) 2487 call pawcprj_alloc(cprj_pos,0,dimcprj) 2488 end if 2489 2490 ! ============================================================================ 2491 ! Loops on positronic bands 2492 2493 do iblock_pos=1,nblock_band_eff_pos 2494 ib_pos=1+(iblock_pos-1)*blocksize 2495 if (any(abs(occ_k_pos(ib_pos:ib_pos+blocksize-1))>tol8)) then 2496 2497 ABI_MALLOC(cwaveg_pos,(2,npw_k_pos*blocksize)) 2498 ABI_MALLOC(cwaveaug_pos,(2,n4,n5,n6*bandpp)) 2499 ABI_MALLOC(denpot_dum,(n4,n5,n6)) 2500 ABI_MALLOC(fofgout_dum,(2,npw_k_pos*blocksize)) 2501 iwavef_pos=(iblock_pos-1)*npw_k_pos*blocksize 2502 cwaveg_pos(:,1:npw_k_pos*blocksize)= & 2503 & cg_k_pos(:,iwavef_pos+1:iwavef_pos+npw_k_pos*blocksize) 2504 2505 ! Get positronic wave function in real space 2506 option=0 2507 if (mpi_enreg%paral_kgb==0) then 2508 weight_pos=occ_k_pos(ib_pos)*wtk_k_pos 2509 call fourwf(1,denpot_dum,cwaveg_pos,fofgout_dum,cwaveaug_pos,& 2510 & gbound_pos,gbound_pos,istwf_k_pos,kg_k_pos,kg_k_pos,& 2511 & dtset%mgfft,mpi_enreg,1,ngfft,npw_k_pos,npw_k_pos,& 2512 & n4,n5,n6,option,tim_fourwf,weight_pos,weight_pos,& 2513 & gpu_option=dtset%gpu_option) 2514 else 2515 call prep_fourwf(denpot_dum,blocksize,cwaveg_pos,cwaveaug_pos,& 2516 & iblock_pos,istwf_k_pos,dtset%mgfft,mpi_enreg,nband_k_pos,& 2517 & bandpp,ngfft,npw_k_pos,n4,n5,n6,occ_k_pos,option,Crystal%ucvol,wtk_k_pos,& 2518 & bandfft_kpt_tab=bandfft_kpt_pos,gpu_option=dtset%gpu_option) 2519 end if 2520 2521 cwaver_pos_block=zero 2522 do ii=1,bandpp 2523 j3=(ii-1)*n3 2524 indx0=1+(ii-1)*cplex*nfft 2525 do i3=1,n3 2526 if (me_fft==fftn3_distrib(i3)) then 2527 indx=indx0+cplex*n1*n2*(ffti3_local(i3)-1) 2528 do i2=1,n2 2529 do i1=1,n1 2530 cwaver_pos_block(indx )=cwaveaug_pos(1,i1,i2,i3+j3) 2531 cwaver_pos_block(indx+1)=cwaveaug_pos(2,i1,i2,i3+j3) 2532 indx=indx+2 2533 end do 2534 end do 2535 end if 2536 end do 2537 end do 2538 ABI_FREE(fofgout_dum) 2539 ABI_FREE(denpot_dum) 2540 ABI_FREE(cwaveaug_pos) 2541 ABI_FREE(cwaveg_pos) 2542 2543 ! At this stage, each band proc has bandpp bands in real space 2544 ! (distributed on FFT procs) 2545 2546 ! ======================================================================== 2547 ! Compute core contribution for this positronic band (PAW only) 2548 2549 if (dtset%usepaw==1) then 2550 do ibpp_pos=1,bandpp 2551 ib_cprj_pos=(iblock_pos-1)*bandpp+ibpp_pos 2552 weight_pos=occ_k_pos(ib_pos+ibpp_pos-1+me_band*bandpp)*wtk_k_pos 2553 ! Calculate the annihilation rate for each core state for state dependent scheme 2554 iatm=0 2555 do itypat=1,dtset%ntypat 2556 mesh_size = pawtab(itypat)%mesh_size 2557 do iat=1,Crystal%nattyp(itypat) 2558 iatm=iatm+1;iatom=Crystal%atindx1(iatm) 2559 ABI_MALLOC(gammastate_c(iatom)%value,(lmncmax(itypat))) 2560 do jlmn=1,lmncmax(itypat) 2561 jln = indlmncor(itypat)%value(5,jlmn) 2562 contrib(:)=zero 2563 ABI_MALLOC(rhocorej,(mesh_size)) 2564 rhocorej(1:mesh_size)=2*phicor(itypat)%value(1:mesh_size,jln)**2 2565 call posratecore(dtset,electronpositron,iatom,dtset%natom,mesh_size,mpi_enreg_seq,& 2566 & 1,pawang,pawrad,pawrhoij_all,pawrhoij_ep_all,pawtab,ratec,rhocorej) 2567 2568 call posratecore(dtset,electronpositron,iatom,dtset%natom,mesh_size,mpi_enreg_seq,& 2569 & 2,pawang,pawrad,pawrhoij_all,pawrhoij_ep_all,pawtab,ratec_ipm,rhocorej) 2570 2571 gammastate_c(iatom)%value(jlmn)=ratec/ratec_ipm 2572 ABI_FREE(rhocorej) 2573 end do 2574 end do 2575 end do 2576 jkpt=0 2577 do ikpt=1,nkpt 2578 if (my_gridtab(ikpt)==0) cycle 2579 jkpt=jkpt+1 2580 do i3=1,n3 2581 ig3=i3-(i3/id3)*n3-1 2582 do i2=1,n2 2583 if (me_fft==fftn2_distrib(i2)) then 2584 j2=ffti2_local(i2) 2585 ig2=i2-(i2/id2)*n2-1 2586 indx=n1*(my_n2*(i3-1)+(j2-1)) 2587 do i1=1,n1 2588 ig1=i1-(i1/id1)*n1-1 2589 indx=indx+1 2590 2591 ! Loop on atoms (type sorted) 2592 iatm=0 2593 do itypat=1,dtset%ntypat 2594 lmn_size = pawtab(itypat)%lmn_size 2595 2596 do iat=1,Crystal%nattyp(itypat) 2597 iatm=iatm+1;iatom=Crystal%atindx1(iatm) 2598 2599 pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+& 2600 & Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+& 2601 & Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt)) 2602 pnorm=dsqrt(dot_product(pcart,pcart)) 2603 pr=dot_product(pcart,Crystal%xcart(:,iatom)) 2604 expipr(1)= cos(two_pi*pr) 2605 expipr(2)=-sin(two_pi*pr) 2606 2607 ! Loop on ij states 2608 do jlmn = 1,lmncmax(itypat) 2609 contrib(:)=zero 2610 do ilmn = 1,lmn_size 2611 radsumnfftc(1)=expipr(1)*radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt)& 2612 & -expipr(2)*radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt) 2613 radsumnfftc(2)=expipr(1)*radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt)& 2614 & +expipr(2)*radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt) 2615 cp_pos(:)=cprj_k_pos(iatom,ib_cprj_pos)%cp(:,ilmn) 2616 contrib(1)=contrib(1)+four_pi*(cp_pos(1)*radsumnfftc(1) & 2617 & -cp_pos(2)*radsumnfftc(2)) 2618 contrib(2)=contrib(2)+four_pi*(cp_pos(1)*radsumnfftc(2) & 2619 & +cp_pos(2)*radsumnfftc(1)) 2620 end do ! end loop over ilmn 2621 ! 2 - electron state weight for 2 spins 2622 rho_moment_core(indx,jkpt) = rho_moment_core(indx,jkpt) & 2623 & +gammastate_c(iatom)%value(jlmn)*2*weight_pos*(contrib(1)**2+contrib(2)**2) 2624 end do ! end loop over jlmn 2625 2626 end do !end loop over atoms 2627 end do !end loop over atom types 2628 2629 end do ! end loop over i1 2630 end if ! end loop over i2 2631 end do 2632 end do ! end loop over i3 2633 end do ! jkpt 2634 end do ! ibpp_pos 2635 end if 2636 2637 ! We now loop over positronic bands inside a block 2638 ! and select occupied ones 2639 do ibpp_pos=1,blocksize 2640 ib_pos=(iblock_pos-1)*blocksize+ibpp_pos 2641 occ_pos=occ_k_pos(ib_pos) 2642 if (abs(occ_pos)>tol8) then 2643 2644 ! Parallelism: dirty trick (broadcast bands) but there should be few positronic bands (~1) 2645 if (nproc_band>1) then 2646 iproc=(ibpp_pos-1)/bandpp 2647 if (me_band==iproc) then 2648 indx=mod((ibpp_pos-1),bandpp)*cplex*nfft 2649 cwaver_pos(1:cplex*nfft)=cwaver_pos_block(indx+1:indx+cplex*nfft) 2650 end if 2651 call xmpi_bcast(cwaver_pos,iproc,mpi_enreg%comm_band,ierr) 2652 if (dtset%usepaw==1) then 2653 if (me_band==iproc) then 2654 indx=mod((ibpp_pos-1),bandpp)*my_nspinor 2655 call pawcprj_copy(cprj_k_pos(:,indx+1:indx+my_nspinor),cprj_pos) 2656 end if 2657 call pawcprj_bcast(cprj_pos,dtset%natom,my_nspinor,dimcprj,0,iproc,& 2658 & mpi_enreg%comm_band,ierr) 2659 end if 2660 else 2661 cwaver_pos(1:cplex*nfft)=cwaver_pos_block(1:cplex*nfft) 2662 if (dtset%usepaw==1) then 2663 call pawcprj_copy(cprj_k_pos(:,(ib_pos-1)*my_nspinor+1:ib_pos*my_nspinor),cprj_pos) 2664 end if 2665 end if 2666 2667 ! ======================================================================== 2668 ! ================ Loop over electronic states =========================== 2669 2670 ! Loop over spins 2671 ibg=0;icg=0;ikg=0;bdtot_index=0 2672 do isppol=1,dtset%nsppol 2673 ! Loop over k points 2674 ikg=0;jkpt=0 2675 do ikpt=1,nkpt 2676 2677 ! Extract data for this kpt_pos 2678 npw_k=npwarr(ikpt) 2679 wtk_k=dtset%wtk(ikpt) 2680 istwf_k=dtset%istwfk(ikpt) 2681 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 2682 nband_cprj_k=nband_k/nproc_band 2683 mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,& 2684 & isppol,mpi_enreg%me_kpt)) 2685 2686 ! Select k-points for current proc 2687 if (mykpt) then 2688 2689 ! Retrieve additional data for this kpt_pos 2690 jkpt=jkpt+1 2691 ABI_MALLOC(occ_k,(nband_k)) 2692 occ_k(:)=occ_ptr(1+bdtot_index:nband_k+bdtot_index) 2693 2694 mcprj_k=0 2695 if (dtset%usepaw==1) then 2696 mcprj_k=my_nspinor*nband_cprj_k 2697 ABI_MALLOC(cprj_k,(dtset%natom,mcprj_k)) 2698 call pawcprj_alloc(cprj_k,0,dimcprj) 2699 call pawcprj_get(Crystal%atindx1,cprj_k,cprj_ptr,dtset%natom,1,ibg,ikpt,iorder_cprj,& 2700 & isppol,mband_cprj,dtset%mkmem,dtset%natom,nband_cprj_k,nband_cprj_k,my_nspinor,& 2701 & dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 2702 end if 2703 2704 if (mpi_enreg%paral_kgb==0) then 2705 ABI_MALLOC(gbound,(2*dtset%mgfft+8,2)) 2706 ABI_MALLOC(kg_k,(3,npw_k)) 2707 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 2708 call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k) 2709 else 2710 jj=mpi_enreg%my_kpttab(ikpt) 2711 bandfft_kpt_el => bandfft_kpt(jj) 2712 end if 2713 2714 ABI_MALLOC(cwaver,(cplex*nfft*bandpp)) 2715 2716 ! ================================================================== 2717 ! Loops on electronic bands 2718 2719 do iblock=1,nblock_band 2720 ib=1+(iblock-1)*blocksize 2721 2722 if (any(abs(occ_k(ib:ib+blocksize-1))>tol8)) then 2723 2724 ! Retrieve electronic wave function 2725 ABI_MALLOC(cwaveg,(2,npw_k*blocksize)) 2726 ABI_MALLOC(cwaveaug,(2,n4,n5,n6*bandpp)) 2727 ABI_MALLOC(denpot_dum,(n4,n5,n6)) 2728 ABI_MALLOC(fofgout_dum,(2,npw_k*blocksize)) 2729 iwavef=(iblock-1)*npw_k*blocksize 2730 cwaveg(:,1:npw_k*blocksize)= & 2731 & cg_ptr(:,icg+iwavef+1:icg+iwavef+npw_k*blocksize) 2732 2733 ! Get electronic wave function in real space 2734 option=0 2735 if (mpi_enreg%paral_kgb==0) then 2736 weight=occ_k(ib)*wtk_k 2737 call fourwf(1,denpot_dum,cwaveg,fofgout_dum,cwaveaug,& 2738 & gbound,gbound,istwf_k,kg_k,kg_k,& 2739 & dtset%mgfft,mpi_enreg,1,ngfft,npw_k,npw_k,& 2740 & n4,n5,n6,option,tim_fourwf,weight,weight,& 2741 & gpu_option=dtset%gpu_option) 2742 else 2743 call prep_fourwf(denpot_dum,blocksize,cwaveg,cwaveaug,& 2744 & iblock,istwf_k,dtset%mgfft,mpi_enreg,nband_k,& 2745 & bandpp,ngfft,npw_k,n4,n5,n6,occ_k,option,Crystal%ucvol,wtk_k,& 2746 & bandfft_kpt_tab=bandfft_kpt_el,gpu_option=dtset%gpu_option) 2747 end if 2748 2749 cwaver=zero 2750 do ii=1,bandpp 2751 j3=(ii-1)*n3 2752 indx0=1+(ii-1)*cplex*nfft 2753 do i3=1,n3 2754 if (me_fft==fftn3_distrib(i3)) then 2755 indx=indx0+cplex*n1*n2*(ffti3_local(i3)-1) 2756 do i2=1,n2 2757 do i1=1,n1 2758 cwaver(indx )=cwaveaug(1,i1,i2,i3+j3) 2759 cwaver(indx+1)=cwaveaug(2,i1,i2,i3+j3) 2760 indx=indx+2 2761 end do 2762 end do 2763 end if 2764 end do 2765 end do 2766 ABI_FREE(fofgout_dum) 2767 ABI_FREE(denpot_dum) 2768 ABI_FREE(cwaveaug) 2769 ABI_FREE(cwaveg) 2770 ! At this stage, each band proc has bandpp bands in real space 2771 ! (distributed on FFT procs) 2772 2773 ! We now loop on the bandpp bands 2774 ! and select occupied ones 2775 do ibpp=1,bandpp 2776 occ_el=occ_k(ib+ibpp-1+me_band*bandpp) 2777 if (abs(occ_el)>tol8) then 2778 2779 ! ============================================================== 2780 ! Compute state-dependent annihilation rate 2781 ! Avoid parallelism over kpt/bands/atoms 2782 gammastate=one;rate_paw=one 2783 if (state_dependent) then 2784 weight=occ_el*wtk_k 2785 ib_cprj=(iblock-1)*bandpp+ibpp 2786 indx=1+(ibpp-1)*cplex*nfft 2787 do ii=1,nfft 2788 rhor_dop_el(ii)=weight*(cwaver(indx)*cwaver(indx)+cwaver(indx+1)*cwaver(indx+1)) 2789 indx=indx+2 2790 end do 2791 if (dtset%usepaw==1) then 2792 do iatom=1,dtset%natom 2793 pawrhoij_dop_el(iatom)%rhoij_=zero 2794 end do 2795 cplex_rhoij=2;if (istwf_k>1) cplex_rhoij=1 2796 use_timerev=(dtset%kptopt>0.and.dtset%kptopt<3) 2797 use_zeromag=(pawrhoij_dop_el(1)%nspden==4.and.dtset%nspden==1) 2798 call pawaccrhoij(Crystal%atindx,cplex_rhoij,cprj_k(:,ib_cprj),& 2799 & cprj_k(:,ib_cprj),0,isppol,dtset%natom,dtset%natom,dtset%nspinor,& 2800 & occ_el,1,pawrhoij_dop_el,use_timerev,use_zeromag,wtk_k) 2801 ! Is it correct to apply symetries here (on a single band)? 2802 ! If not, call pawrhoij_symrhoij with nsym=1 2803 call pawrhoij_symrhoij(pawrhoij_dop_el,pawrhoij_dop_el,1,Crystal%gprimd,& 2804 & Crystal%indsym,0,dtset%natom,Crystal%nsym,dtset%ntypat,1,pawang,-10001,& 2805 & pawtab,Crystal%rprimd,Crystal%symafm,Crystal%symrec,dtset%typat) 2806 end if 2807 ! Has to call poslifetime in sequential because we are in a parallel section 2808 ! Only FFT parallelism is allowed 2809 call poslifetime(dtset,electronpositron,Crystal%gprimd,dtset%natom,mpi_enreg_seq,n3xccc,& 2810 & nfft,ngfft,nhat,2,pawang,pawrad,pawrhoij_all,pawtab,rate,rate_paw,rhor,Crystal%ucvol,xccc3d,& 2811 & rhor_dop_el=rhor_dop_el,pawrhoij_dop_el=pawrhoij_dop_el,pawrhoij_ep=pawrhoij_ep_all) 2812 call poslifetime(dtset,electronpositron,Crystal%gprimd,dtset%natom,mpi_enreg_seq,n3xccc,& 2813 & nfft,ngfft,nhat,3,pawang,pawrad,pawrhoij_all,pawtab,rate_ipm,rate_paw_ipm,rhor,Crystal%ucvol,xccc3d,& 2814 & rhor_dop_el=rhor_dop_el,pawrhoij_dop_el=pawrhoij_dop_el,pawrhoij_ep=pawrhoij_ep_all) 2815 gammastate=rate/rate_ipm 2816 rate_paw=rate_paw/rate_paw_ipm 2817 end if 2818 2819 ! ============================================================== 2820 ! Compute plane-wave contribution to momentum distribution 2821 2822 ! Compute Psi^+(r) * Psi^-(r) * gamma(r) in real space 2823 rho_contrib(:)=zero 2824 indx=(ibpp-1)*cplex*nfft 2825 if (cplex==2) then 2826 do jj=1,nfft 2827 ii=2*jj-1 2828 rho_contrib(ii) =sqrt(gamma(jj,2))*(cwaver_pos(ii)*cwaver(indx+ii)& 2829 & -wf_fact*cwaver_pos(ii+1)*cwaver(indx+ii+1)) 2830 rho_contrib(ii+1)=sqrt(gamma(jj,2))*(cwaver_pos(ii)*cwaver(indx+ii+1) & 2831 & +wf_fact*cwaver_pos(ii+1)*cwaver(indx+ii)) 2832 end do 2833 else 2834 do ii=1,nfft 2835 rho_contrib(ii)=sqrt(gamma(ii,2))*cwaver_pos(ii)*cwaver(indx+ii) 2836 end do 2837 end if 2838 2839 ! FFT of (Psi+.Psi-.gamma) to get Intg[(Psi+.Psi-.gamma).exp(-igr)] 2840 call fourdp(cplex,rho_contrib_g,rho_contrib,-1,mpi_enreg,nfft,1,ngfft,& 2841 & tim_fourdp) 2842 2843 rho_pw(1:nfft,jkpt)=rho_pw(1:nfft,jkpt) +gammastate*occ_el*occ_pos & 2844 & *(rho_contrib_g(1,1:nfft)**2+rho_contrib_g(2,1:nfft)**2) 2845 2846 ! ============================================================== 2847 ! Compute PAW on-site contribution to momentum distribution 2848 2849 if (dtset%usepaw==1) then 2850 2851 rho_contrib_paw1(:,:)= zero 2852 rho_contrib_paw2(:,:)= zero 2853 rho_contrib_paw3(:,:)= zero 2854 2855 ib_cprj=(iblock-1)*bandpp+ibpp 2856 2857 ! Loop on moments 2858 indx=0 2859 do i3=1,n3 2860 ig3=i3-(i3/id3)*n3-1 2861 do i2=1,n2 2862 if (me_fft==fftn2_distrib(i2)) then 2863 j2=ffti2_local(i2) 2864 ig2=i2-(i2/id2)*n2-1 2865 indx=n1*(my_n2*(i3-1)+(j2-1)) 2866 do i1=1,n1 2867 ig1=i1-(i1/id1)*n1-1 2868 indx=indx+1 2869 2870 pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+& 2871 & Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+& 2872 & Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt)) 2873 pnorm=dsqrt(dot_product(pcart,pcart)) 2874 2875 ! Loop on atoms (type-sorted) 2876 iatm=0 2877 do itypat=1,dtset%ntypat 2878 lmn_size=pawtab(itypat)%lmn_size 2879 lmn2_size=pawtab(itypat)%lmn2_size 2880 ABI_MALLOC(radsumnfft1,(2,lmn2_size)) 2881 ABI_MALLOC(radsumnfft2,(2,lmn2_size)) 2882 ABI_MALLOC(radsumnfft3,(2,lmn2_size)) 2883 2884 do iat=1,Crystal%nattyp(itypat) 2885 iatm=iatm+1;iatom=Crystal%atindx1(iatm) 2886 2887 pr=dot_product(pcart,Crystal%xcart(:,iatom)) 2888 expipr(1)= cos(two_pi*pr) 2889 expipr(2)=-sin(two_pi*pr) 2890 2891 do klmn=1,lmn2_size 2892 radsumnfft1(1,klmn)=expipr(1)*radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt)& 2893 & -expipr(2)*radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt) 2894 radsumnfft1(2,klmn)=expipr(1)*radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt)& 2895 & +expipr(2)*radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt) 2896 radsumnfft2(1,klmn)=expipr(1)*radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt)& 2897 & -expipr(2)*radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt) 2898 radsumnfft2(2,klmn)=expipr(1)*radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt)& 2899 & +expipr(2)*radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt) 2900 radsumnfft3(1,klmn)=expipr(1)*radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt)& 2901 & -expipr(2)*radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt) 2902 radsumnfft3(2,klmn)=expipr(1)*radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt)& 2903 & +expipr(2)*radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt) 2904 end do 2905 2906 ! Loop on ij states 2907 do ilmn = 1, lmn_size 2908 i0lmn = ilmn*(ilmn-1)/2 2909 do jlmn = 1, lmn_size 2910 klmn = i0lmn+jlmn 2911 if (jlmn>ilmn) then 2912 i0lmn=jlmn*(jlmn-1)/2; klmn=i0lmn+ilmn 2913 end if 2914 ! Transform 3-dimentional radsum to 1-dimentional radsumnfft 2915 cp(:)=cprj_k(iatom,ib_cprj)%cp(:,ilmn) 2916 cp_pos(:)=cprj_pos(iatom,1)%cp(:,jlmn) 2917 cp11= cp(1)*cp_pos(1) 2918 cp22= cp(2)*cp_pos(2)*wf_fact 2919 cp12= cp(1)*cp_pos(2)*wf_fact 2920 cp21= cp(2)*cp_pos(1) 2921 cpr=cp11-cp22 ; cpi=cp12+cp21 2922 rho_contrib_paw1(1,indx) = rho_contrib_paw1(1,indx) & 2923 & + four_pi*(cpr*radsumnfft1(1,klmn)-cpi*radsumnfft1(2,klmn)) 2924 rho_contrib_paw1(2,indx) = rho_contrib_paw1(2,indx) & 2925 & + four_pi*(cpr*radsumnfft1(2,klmn)+cpi*radsumnfft1(1,klmn)) 2926 rho_contrib_paw2(1,indx) = rho_contrib_paw2(1,indx) & 2927 & + four_pi*(cpr*radsumnfft2(1,klmn)-cpi*radsumnfft2(2,klmn)) 2928 rho_contrib_paw2(2,indx) = rho_contrib_paw2(2,indx) & 2929 & + four_pi*(cpr*radsumnfft2(2,klmn)+cpi*radsumnfft2(1,klmn)) 2930 rho_contrib_paw3(1,indx) = rho_contrib_paw3(1,indx) & 2931 & + four_pi*(cpr*radsumnfft3(1,klmn)-cpi*radsumnfft3(2,klmn)) 2932 rho_contrib_paw3(2,indx) = rho_contrib_paw3(2,indx) & 2933 & + four_pi*(cpr*radsumnfft3(2,klmn)+cpi*radsumnfft3(1,klmn)) 2934 end do ! end loop over jlmn 2935 end do ! end loop over ilmn 2936 2937 end do !end loop over atoms 2938 2939 ABI_FREE(radsumnfft1) 2940 ABI_FREE(radsumnfft2) 2941 ABI_FREE(radsumnfft3) 2942 end do !end loop over atom types 2943 2944 rho_moment_v1(indx,jkpt) = rho_moment_v1(indx,jkpt) & 2945 & +occ_el*occ_pos & 2946 & *(gammastate*(rho_contrib_g(1,indx)**2+rho_contrib_g(2,indx)**2) & 2947 & +rate_paw*(rho_contrib_paw1(1,indx)**2+rho_contrib_paw1(2,indx)**2 & 2948 & -rho_contrib_paw2(1,indx)**2-rho_contrib_paw2(2,indx)**2)) 2949 rho_moment_v2(indx,jkpt) = rho_moment_v2(indx,jkpt) & 2950 & +occ_el*occ_pos*gammastate & 2951 & *((rho_contrib_g(1,indx)+rho_contrib_paw3(1,indx))**2+& 2952 & (rho_contrib_g(2,indx)+rho_contrib_paw3(2,indx))**2) 2953 2954 end do ! end loop over i1 2955 2956 end if ! end loop over i2 2957 end do 2958 end do ! end loop over i3 2959 2960 end if ! PAW 2961 2962 ! ================================================================ 2963 ! End loops on electronic bands 2964 2965 end if ! occ>1.e-8 2966 end do ! ibpp 2967 end if ! occ_block>1.e-8 2968 end do ! iblock 2969 2970 ! End loops over k points and spins (electrons) 2971 icg = icg + npw_k*my_nspinor*nband_k 2972 ibg = ibg + my_nspinor*nband_cprj_k 2973 ikg = ikg + npw_k 2974 2975 ABI_FREE(cwaver) 2976 ABI_FREE(occ_k) 2977 if (mpi_enreg%paral_kgb==0) then 2978 ABI_FREE(kg_k) 2979 ABI_FREE(gbound) 2980 else 2981 nullify(bandfft_kpt_el) 2982 end if 2983 if (dtset%usepaw==1) then 2984 call pawcprj_free(cprj_k) 2985 ABI_FREE(cprj_k) 2986 end if 2987 2988 end if ! mykpt 2989 bdtot_index=bdtot_index+nband_k 2990 end do ! ikpt 2991 end do ! isppol 2992 2993 ! ================================================================ 2994 ! End loops on positronic bands 2995 2996 end if ! occ>1.e-8 2997 end do ! ibpp_pos 2998 end if ! occ(block)>1.e-8 2999 end do ! iblock_pos 3000 3001 ! End loop over k points (positron) 3002 if (mykpt_pos) then 3003 icg_pos = icg_pos + npw_k_pos*my_nspinor*nband_k_pos 3004 ibg_pos = ibg_pos + my_nspinor*nband_cprj_k_pos 3005 ikg_pos = ikg_pos + npw_k_pos 3006 end if 3007 bdtot_index_pos=bdtot_index_pos+nband_k_pos 3008 3009 ABI_FREE(cwaver_pos) 3010 ABI_FREE(cwaver_pos_block) 3011 ABI_FREE(cg_k_pos) 3012 ABI_FREE(occ_k_pos) 3013 if (mpi_enreg%paral_kgb==0) then 3014 ABI_FREE(kg_k_pos) 3015 ABI_FREE(gbound_pos) 3016 else if (mykpt_pos) then 3017 nullify(bandfft_kpt_pos) 3018 else 3019 call bandfft_kpt_destroy(bandfft_kpt_pos) 3020 ABI_FREE(bandfft_kpt_pos) 3021 end if 3022 if (dtset%usepaw==1) then 3023 call pawcprj_free(cprj_pos) 3024 ABI_FREE(cprj_pos) 3025 call pawcprj_free(cprj_k_pos) 3026 ABI_FREE(cprj_k_pos) 3027 end if 3028 3029 end do ! ikpt_pos 3030 3031 !================================================================ 3032 !Final computations and printing 3033 3034 !In case of parallelism, sum over the communicator(s) 3035 if (nproc_band>1) then 3036 ABI_MALLOC(mpibuf,(3*nfft,my_ngrid)) 3037 do jkpt=1,my_ngrid 3038 mpibuf( 1: nfft,jkpt)=rho_moment_v1(1:nfft,jkpt) 3039 mpibuf( nfft+1:2*nfft,jkpt)=rho_moment_v2(1:nfft,jkpt) 3040 mpibuf(2*nfft+1:3*nfft,jkpt)=rho_pw (1:nfft,jkpt) 3041 end do 3042 call xmpi_sum(mpibuf,mpi_enreg%comm_band,ierr) 3043 do jkpt=1,my_ngrid 3044 rho_moment_v1(1:nfft,jkpt)=mpibuf( 1: nfft,jkpt) 3045 rho_moment_v2(1:nfft,jkpt)=mpibuf( nfft+1:2*nfft,jkpt) 3046 rho_pw(1:nfft,jkpt) =mpibuf(2*nfft+1:3*nfft,jkpt) 3047 end do 3048 ABI_FREE(mpibuf) 3049 end if 3050 if (dtset%usepaw==1) then 3051 call xmpi_sum(rho_moment_core,mpi_enreg%comm_band,ierr) 3052 end if 3053 3054 !Add valence and core contributions 3055 if (dtset%usepaw==1) then 3056 if (dtset%nsppol==2.and.my_nsppol==1) rho_moment_core(:,:)=half*rho_moment_core(:,:) 3057 rho_moment_v1(:,:)=rho_moment_v1(:,:)+rho_moment_core(:,:) 3058 rho_moment_v2(:,:)=rho_moment_v2(:,:)+rho_moment_core(:,:) 3059 end if 3060 3061 units_=pi*(one/InvFineStruct)**3/Time_Sec/1.e12_dp/electronpositron%posocc 3062 scale_=(two_pi**2)/(Crystal%ucvol**two_thirds) 3063 3064 !Integrate rho_moment over p 3065 buf(1)=sum(rho_moment_v1(1:nfft,1:my_ngrid)) 3066 buf(2)=sum(rho_moment_v2(1:nfft,1:my_ngrid)) 3067 buf(3)=sum(rho_moment_core(1:nfft,1:my_ngrid)) 3068 buf(4)=sum(rho_pw(1:nfft,1:my_ngrid)) 3069 call xmpi_sum(buf,mpi_enreg%comm_kpt,ierr) 3070 call xmpi_sum(buf,mpi_enreg%comm_fft,ierr) 3071 lambda_v1=buf(1)*units_/Crystal%ucvol/nkpt 3072 lambda_v2=buf(2)*units_/Crystal%ucvol/nkpt 3073 lambda_core=buf(3)*units_/Crystal%ucvol/nkpt 3074 lambda_pw=buf(4)*units_/Crystal%ucvol/nkpt 3075 3076 !Write result in _DOPPLER file 3077 !Requires MPI-IO if nproc_fft>1 3078 if (me_band==0) then 3079 if (me_kpt==0) then 3080 filename_dop=trim(dtfil%filnam_ds(4))//'_DOPPLER' 3081 vec=sqrt(dot_product(Crystal%gprimd(:,3),Crystal%gprimd(:,3))) 3082 ABI_MALLOC(pcart_k,(3,nfft)) 3083 ABI_MALLOC(rho_moment_k,(nfft)) 3084 if (dtset%nsppol==2) then 3085 ABI_MALLOC(rho_moment_k2,(nfft)) 3086 end if 3087 if (accessfil==IO_MODE_FORTRAN) then ! >>>>> Fortran access 3088 ! Open file and write first line 3089 ierr=open_file(filename_dop,msg,newunit=unit_doppler,form='unformatted') 3090 write(unit_doppler) nfft,nkpt,Crystal%ucvol,Crystal%rprimd(:,:) 3091 else ! >>>>> MPI-IO access 3092 unit_doppler=get_unit() 3093 ! Open file and write first line 3094 call WffOpen(IO_MODE_MPI,mpi_enreg%comm_fft,filename_dop,ierr,wff,0,me_fft,unit_doppler) 3095 if (me_fft==0) then 3096 call xderiveWRecInit(wff,ierr) 3097 call xderiveWrite(wff,n1*n2*n3,ierr) 3098 call xderiveWrite(wff,nkpt,ierr) 3099 call xderiveWrite(wff,Crystal%ucvol,ierr) 3100 call xderiveWrite(wff,Crystal%rprimd(:,:),ierr) 3101 call xderiveWRecEnd(wff,ierr) 3102 else 3103 call xmoveOff(wff,n_int=2,n_dp=10,n_mark=2) 3104 end if 3105 ! Store table of FFT points treated by current proc 3106 ABI_MALLOC(my_ffttab,(nfft)) 3107 my_ffttab=0 3108 do i3=1,n3 3109 do i2=1,n2 3110 if (me_fft==fftn2_distrib(i2)) then 3111 indx0=n1*(n2*(i3-1)+(i2-1)) 3112 indx=n1*(my_n2*(i3-1)+(ffti2_local(i2)-1)) 3113 my_ffttab(indx+1:indx+n1)=(/(indx0+ii,ii=1,n1)/) 3114 end if 3115 end do 3116 end do 3117 ABI_MALLOC(mpibuf,(1,nfft)) 3118 end if 3119 end if 3120 3121 jkpt=0 3122 do ikpt=1,nkpt 3123 if (nproc_spkpt==1) then 3124 rho_moment_k(1:nfft)=rho_moment_v2(1:nfft,ikpt) 3125 else 3126 if (my_gridtab(ikpt)/=0) jkpt=jkpt+1 3127 if (me_kpt==0) then 3128 if (my_gridtab(ikpt)==0) then 3129 tag=ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,1) 3130 call xmpi_recv(rho_moment_k,iproc,tag,mpi_enreg%comm_kpt,ierr) 3131 if (dtset%nsppol==2) then 3132 tag=2*ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,2) 3133 call xmpi_recv(rho_moment_k2,iproc,tag,mpi_enreg%comm_kpt,ierr) 3134 rho_moment_k(1:nfft)=rho_moment_k(1:nfft)+rho_moment_k2(1:nfft) 3135 end if 3136 else if (any(mpi_enreg%my_isppoltab(:)==1)) then 3137 rho_moment_k(1:nfft)=rho_moment_v2(1:nfft,jkpt) 3138 if (dtset%nsppol==2) then 3139 ii=2;if (mpi_enreg%my_isppoltab(2)==1) ii=1 3140 tag=ii*ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,ii) 3141 call xmpi_recv(rho_moment_k2,iproc,tag,mpi_enreg%comm_kpt,ierr) 3142 rho_moment_k(1:nfft)=rho_moment_k(1:nfft)+rho_moment_k2(1:nfft) 3143 end if 3144 end if 3145 else if (my_gridtab(ikpt)/=0) then 3146 if (mpi_enreg%my_isppoltab(1)==1) then 3147 tag=ikpt 3148 call xmpi_send(rho_moment_v2(1:nfft,jkpt),0,tag,mpi_enreg%comm_kpt,ierr) 3149 end if 3150 if (dtset%nsppol==2.and.mpi_enreg%my_isppoltab(2)==1) then 3151 tag=2*ikpt 3152 call xmpi_send(rho_moment_v2(1:nfft,jkpt),0,tag,mpi_enreg%comm_kpt,ierr) 3153 end if 3154 end if 3155 end if ! nproc_spkpt>1 3156 if (me_kpt==0) then 3157 indx=0 3158 do i3=1,n3 3159 ig3=i3-(i3/id3)*n3-1 3160 do i2=1,n2 3161 if (me_fft/=fftn2_distrib(i2)) cycle 3162 ig2=i2-(i2/id2)*n2-1 3163 do i1=1,n1 3164 ig1=i1-(i1/id1)*n1-1 3165 indx=indx+1 3166 pcart_k(:,indx)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+& 3167 & Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+& 3168 & Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt)) 3169 end do 3170 end do 3171 end do 3172 if (accessfil==IO_MODE_FORTRAN) then 3173 write(unit_doppler) pcart_k(1:3,1:nfft),rho_moment_k(1:nfft) 3174 else 3175 mpibuf(1,1:nfft)=rho_moment_k(1:nfft) 3176 call xderiveWRecInit(wff,ierr) 3177 call xderiveWrite(wff,pcart_k,3,nfft,mpi_enreg%comm_fft,my_ffttab,ierr) 3178 call xderiveWrite(wff,mpibuf ,1,nfft,mpi_enreg%comm_fft,my_ffttab,ierr) 3179 call xderiveWRecEnd(wff,ierr) 3180 end if 3181 end if 3182 end do 3183 if (me_kpt==0) then 3184 ABI_FREE(pcart_k) 3185 ABI_FREE(rho_moment_k) 3186 if (dtset%nsppol==2) then 3187 ABI_FREE(rho_moment_k2) 3188 end if 3189 if (accessfil==IO_MODE_FORTRAN) then 3190 ierr=close_unit(unit_doppler,msg) 3191 else 3192 call WffClose(wff,ierr) 3193 ABI_FREE(my_ffttab) 3194 ABI_FREE(mpibuf) 3195 end if 3196 end if 3197 end if ! me_band==0 3198 3199 !Write results 3200 write(msg,'(7a)') & 3201 & ' Computation of electron-positron pairs momentum distribution completed.',ch10,& 3202 & '-File ',trim(filename_dop),' has been created.',ch10,& 3203 & '-Use ~abinit/scripts/post_processing/posdopspectra.F90 to process it.' 3204 call wrtout(ab_out,msg,'COLL') 3205 call wrtout(std_out,msg,'COLL') 3206 msg=' Some e-p annihilation rates (ns-1) obtained by integration of e-p pairs momentum distribution:' 3207 call wrtout(std_out,msg,'COLL') 3208 write(msg,'(a,es22.12,3(2a,es22.12))') & 3209 & ' Lambda (from module of sum of PAW contrib.) = ',lambda_v2*1000._dp,ch10,& 3210 & ' = lambda_core: ',lambda_core*1000._dp,ch10,& 3211 & ' +lambda_pw : ',lambda_pw*1000._dp,ch10,& 3212 & ' +lambda_paw : ',(lambda_v2-lambda_core-lambda_pw)*1000._dp 3213 call wrtout(std_out,msg,'COLL') 3214 write(msg,'(4(a,es22.12,a))') & 3215 & ' Lambda (from sum of modules of PAW contrib.) = ',lambda_v1*1000._dp,ch10,& 3216 & ' = lambda_core: ',lambda_core*1000._dp,ch10,& 3217 & ' +lambda_pw : ',lambda_pw*1000._dp,ch10,& 3218 & ' +lambda_paw : ',(lambda_v1-lambda_core-lambda_pw)*1000._dp,ch10 3219 call wrtout(std_out,msg,'COLL') 3220 write(msg,'(4a,es22.12,2a)') ch10,& 3221 & ' Annihilation rate obtained from integration of e-p pairs momentum distribution:',ch10,& 3222 & ' lambda=',lambda_v2*1000._dp,' ns-1',ch10 3223 call wrtout(ab_out,msg,'COLL') 3224 3225 !Deallocate remaining memory 3226 ABI_FREE(my_gridtab) 3227 ABI_FREE(rho_pw) 3228 ABI_FREE(rho_moment_v1) 3229 ABI_FREE(rho_moment_v2) 3230 ABI_FREE(rho_moment_core) 3231 ABI_FREE(rho_contrib) 3232 ABI_FREE(rho_contrib_g) 3233 ABI_FREE(rho_contrib_paw1) 3234 ABI_FREE(rho_contrib_paw2) 3235 ABI_FREE(rho_contrib_paw3) 3236 if (state_dependent) then 3237 call unset_mpi_enreg_fft(mpi_enreg_seq) 3238 call destroy_mpi_enreg(mpi_enreg_seq) 3239 ABI_FREE(rhor_dop_el) 3240 if (dtset%usepaw==1) then 3241 call pawrhoij_free(pawrhoij_dop_el) 3242 ABI_FREE(pawrhoij_dop_el) 3243 if (mpi_enreg%my_natom<dtset%natom) then 3244 call pawrhoij_free(pawrhoij_all) 3245 call pawrhoij_free(pawrhoij_ep_all) 3246 ABI_FREE(pawrhoij_all) 3247 ABI_FREE(pawrhoij_ep_all) 3248 end if 3249 end if 3250 end if 3251 3252 ABI_FREE(gamma) 3253 3254 if (dtset%usepaw==1.and.(.not.include_nhat_in_gamma)) then 3255 ABI_FREE(rhor_) 3256 ABI_FREE(rhor_ep_) 3257 end if 3258 3259 if (dtset%usepaw==1) then 3260 ABI_FREE(nphicor) 3261 ABI_FREE(lmncmax) 3262 do itypat=1,dtset%ntypat 3263 if (allocated(phicor(itypat)%value)) then 3264 ABI_FREE(phicor(itypat)%value) 3265 end if 3266 if (allocated(indlmncor(itypat)%value)) then 3267 ABI_FREE(indlmncor(itypat)%value) 3268 end if 3269 if (allocated(radsumc(itypat)%value)) then 3270 ABI_FREE(radsumc(itypat)%value) 3271 end if 3272 if (allocated(radsum1(itypat)%value)) then 3273 ABI_FREE(radsum1(itypat)%value) 3274 end if 3275 if (allocated(radsum2(itypat)%value)) then 3276 ABI_FREE(radsum2(itypat)%value) 3277 end if 3278 if (allocated(radsum3(itypat)%value)) then 3279 ABI_FREE(radsum3(itypat)%value) 3280 end if 3281 end do 3282 do iatom=1,dtset%natom 3283 if (allocated(gammastate_c(iatom)%value)) then 3284 ABI_FREE(gammastate_c(iatom)%value) 3285 end if 3286 end do 3287 ABI_FREE(phicor) 3288 ABI_FREE(indlmncor) 3289 ABI_FREE(radsumc) 3290 ABI_FREE(radsum1) 3291 ABI_FREE(radsum2) 3292 ABI_FREE(radsum3) 3293 ABI_FREE(gammastate_c) 3294 end if 3295 3296 DBG_EXIT("COLL") 3297 3298 end subroutine posdoppler
ABINIT/poslifetime [ Functions ]
NAME
poslifetime
FUNCTION
Calculate the positron lifetime
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset | nspden=number of spin-density components | ntypat=number of atom types | paral_kgb=flag controlling (k,g,bands) parallelization | pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) | usepaw=flag for PAW gprimd(3,3)= dimensional reciprocal space primitive translations mpi_enreg= information about MPI parallelization my_natom=number of atoms treated by current processor n3xccc= dimension of the xccc3d array (0 or nfft). nfft= number of FFT grid points ngfft(18)= contain all needed information about 3D FFT nhat(nfft,nspden)=charge compensation density (content depends on electronpositron%particle) option= if 1, calculate positron lifetime for whole density if 2, calculate positron lifetime for given state if 3, calculate positron lifetime for given state with IPM pawang <type(pawang)>=paw angular mesh and related data pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data rhor(nfft,nspden)=total electron/positron density (content depends on electronpositron%particle) ucvol=unit cell volume in bohr**3. xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 ===== Optional arguments, used only if option>1 ===== pawrhoij_dop_el(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies of one state rhor_dop_el(nfft)=electron density of given state for the state dependent scheme ===== Optional argument ===== pawrhoij_ep(my_natom*usepaw) <type(pawrhoij_type)>= atomic occupancies to be used in place of electronpositron%pawrhoij_ep
OUTPUT
rate= annihilation rate of a given state needed for state dependent scheme for doppler broadening
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
SOURCE
911 subroutine poslifetime(dtset,electronpositron,gprimd,my_natom,mpi_enreg,n3xccc,nfft,ngfft,nhat,& 912 & option,pawang,pawrad,pawrhoij,pawtab,rate,rate_paw,rhor,ucvol,xccc3d,& 913 & rhor_dop_el,pawrhoij_dop_el,pawrhoij_ep) ! optional arguments 914 915 !Arguments ------------------------------------ 916 !scalars 917 integer,intent(in) :: my_natom,n3xccc,nfft,option 918 real(dp),intent(in) :: ucvol 919 real(dp),intent(out) :: rate,rate_paw 920 type(dataset_type), intent(in) :: dtset 921 type(electronpositron_type),pointer :: electronpositron 922 type(MPI_type),intent(in) :: mpi_enreg 923 type(pawang_type), intent(in) :: pawang 924 !arrays 925 integer,intent(in) :: ngfft(18) 926 real(dp),intent(in) :: gprimd(3,3),nhat(nfft,dtset%nspden*dtset%usepaw),xccc3d(n3xccc) 927 real(dp),intent(in),target :: rhor(nfft,dtset%nspden) 928 real(dp),optional,intent(in) :: rhor_dop_el(nfft) 929 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw) 930 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*dtset%usepaw) 931 type(pawrhoij_type),optional,intent(in) :: pawrhoij_dop_el(my_natom*dtset%usepaw) 932 type(pawrhoij_type),optional,target,intent(in) :: pawrhoij_ep(my_natom*dtset%usepaw) 933 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 934 935 !Local variables------------------------------- 936 !scalars 937 integer :: cplex,iatom,ierr,ifft,igam,ii,ilm,ilm1,ilm2,iloop,ipt,ir,isel 938 integer :: itypat,iwarn,iwarnj,iwarnp,lm_size,lmn2_size,mesh_size 939 integer :: nfftot,ngamma,ngr,ngrad,nspden_ep,opt_dens,usecore 940 logical,parameter :: include_nhat_in_gamma=.false. 941 real(dp),parameter :: delta=1.d-4 942 real(dp) :: fact,fact2,intg 943 real(dp) :: lambda_core ,lambda_core_ipm ,lambda ,lambda_ipm 944 real(dp) :: lambda_core_paw,lambda_core_paw_ipm,lambda_paw,lambda_paw_ipm 945 real(dp) :: lifetime,lifetime_ipm,nbec,nbev,nbp,rdum,sqfpi,units 946 character(len=500) :: msg 947 !arrays 948 integer,allocatable :: igamma(:) 949 logical,allocatable :: lmselect(:),lmselect_ep(:),lmselect_dum(:) 950 real(dp) :: mpibuf(4) 951 real(dp),parameter :: qphon(3)=(/zero,zero,zero/),lsign(2)=(/one,-one/) 952 real(dp),allocatable :: d1gam(:,:,:),d2gam(:,:,:),ff(:),gam_(:,:,:),gamma(:,:),gammam(:,:,:),gg(:,:) 953 real(dp),allocatable :: grhocore2(:),grhocor2_(:),grhoe2(:),grho2_(:) 954 real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:),nhat1_j(:,:,:) 955 real(dp),allocatable :: rho_(:),rho_ep_(:),rho1(:,:,:),rho1_ep(:,:,:),rho1_j(:,:,:) 956 real(dp),allocatable :: rhoarr1(:),rhoarr1_ep(:),rhoarr1_j(:),rhoarr2(:) 957 real(dp),allocatable :: rhocore(:),rhocor_(:),rhoe(:,:),rhop(:,:),rhor_dop_el_(:) 958 real(dp),allocatable :: rhosph(:),rhosph_ep(:),rhosph_j(:),rhotot(:,:),rhotot_ep(:,:) 959 real(dp),allocatable :: rhotot_j(:,:),trho1(:,:,:),trho1_ep(:,:,:),trho1_j(:,:,:) 960 real(dp),allocatable :: v1sum(:,:),v2sum(:,:,:) 961 real(dp),pointer :: rhor_(:,:),rhor_ep_(:,:) 962 type(pawrhoij_type),pointer :: pawrhoij_ep_(:) 963 964 ! ************************************************************************* 965 966 DBG_ENTER("COLL") 967 968 !Tests for developers 969 if (.not.associated(electronpositron)) then 970 msg='electronpositron variable must be associated!' 971 ABI_BUG(msg) 972 end if 973 if (option/=1) then 974 if ((.not.present(rhor_dop_el)).or.(.not.present(pawrhoij_dop_el))) then 975 msg='when option/=1, rhor_dop_el and pawrhoij_dop_el must be present!' 976 ABI_BUG(msg) 977 end if 978 end if 979 980 ! This to avoid using unitialized variables. 981 lambda_core = zero; lambda_paw = zero; lambda_core_paw = zero 982 983 !Constants 984 fact=0.0 985 cplex=1;nspden_ep=1 986 usecore=n3xccc/nfft 987 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 988 ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2 989 iwarn=0;iwarnj=0;iwarnp=1 990 sqfpi=sqrt(four_pi) 991 992 !Compatibility tests 993 if (electronpositron%particle==EP_NOTHING) then 994 msg='Not valid for electronpositron%particle=NOTHING!' 995 ABI_BUG(msg) 996 end if 997 if (electronpositron%nfft/=nfft) then 998 msg='nfft/=electronpositron%nfft!' 999 ABI_BUG(msg) 1000 end if 1001 if (dtset%usepaw==1) then 1002 if(dtset%pawxcdev==0.and.ngrad==2) then 1003 msg='GGA is not implemented for pawxcdev=0 (use dtset%pawxcdev/=0)!' 1004 ABI_BUG(msg) 1005 end if 1006 end if 1007 1008 !Select type(s) of enhancement factor 1009 if ((electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3).and.option==1) then 1010 ngamma=2 1011 ABI_MALLOC(igamma,(ngamma)) 1012 igamma(1)=1;igamma(2)=2 1013 else 1014 ngamma=1 1015 ABI_MALLOC(igamma,(ngamma)) 1016 if (electronpositron%ixcpositron==-1) igamma(1)=0 1017 if (electronpositron%ixcpositron== 2) igamma(1)=4 1018 if (electronpositron%ixcpositron==11.or.electronpositron%ixcpositron==31) igamma(1)=3 1019 if (electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3) igamma(1)=2 1020 end if 1021 1022 !Select density according to nhat choice 1023 if (dtset%usepaw==0.or.include_nhat_in_gamma) then 1024 rhor_ => rhor 1025 rhor_ep_ => electronpositron%rhor_ep 1026 else 1027 ABI_MALLOC(rhor_,(nfft,dtset%nspden)) 1028 ABI_MALLOC(rhor_ep_,(nfft,dtset%nspden)) 1029 rhor_=rhor-nhat 1030 rhor_ep_=electronpositron%rhor_ep-electronpositron%nhat_ep 1031 end if 1032 1033 !Eventually overwrite electronpositron%pawrhoij_ep 1034 if (present(pawrhoij_ep)) then 1035 pawrhoij_ep_ => pawrhoij_ep 1036 else 1037 pawrhoij_ep_ => electronpositron%pawrhoij_ep 1038 end if 1039 1040 !Loop on different enhancement factors 1041 do igam=1,ngamma 1042 1043 ! Compute electron-positron annihilation rate using pseudo densities (plane waves) 1044 ! ---------------------------------------------------------------------------------------- 1045 1046 ! Select the densities and make them positive 1047 ABI_MALLOC(rhoe,(nfft,nspden_ep)) 1048 ABI_MALLOC(rhop,(nfft,nspden_ep)) 1049 if (electronpositron%particle==EP_ELECTRON) then 1050 rhoe(:,1)=rhor_ep_(:,1);rhop(:,1)=rhor_(:,1) 1051 else if (electronpositron%particle==EP_POSITRON) then 1052 rhoe(:,1)=rhor_(:,1);rhop(:,1)=rhor_ep_(:,1) 1053 end if 1054 call mkdenpos(iwarn ,nfft,nspden_ep,1,rhoe,dtset%xc_denpos) 1055 call mkdenpos(iwarnp,nfft,nspden_ep,1,rhop,dtset%xc_denpos) 1056 if (option/=1) then 1057 ABI_MALLOC(rhor_dop_el_,(nfft)) 1058 rhor_dop_el_(:)=rhor_dop_el(:) 1059 call mkdenpos(iwarnp,nfft,1,1,rhor_dop_el_,dtset%xc_denpos) 1060 end if 1061 1062 ! Compute enhancement factor at each FFT grid point 1063 ! gamma(:,1): using total electronic density 1064 ! gamma(:,2): using valence electronic density 1065 ABI_MALLOC(gamma,(nfft,2)) 1066 if (option==1.or.option==2) then 1067 call gammapositron_fft(electronpositron,gamma,gprimd,igamma(igam),mpi_enreg,& 1068 & n3xccc,nfft,ngfft,rhoe,rhop,xccc3d) 1069 else 1070 gamma=one 1071 end if 1072 1073 ! Compute positron annihilation rates 1074 lambda =zero;lambda_ipm =zero 1075 lambda_core=zero;lambda_core_ipm=zero 1076 if (option==1) then 1077 do ifft=1,nfft 1078 lambda =lambda +rhop(ifft,1)*rhoe(ifft,1)*gamma(ifft,1) 1079 lambda_ipm=lambda_ipm+rhop(ifft,1)*rhoe(ifft,1)*gamma(ifft,2) 1080 end do 1081 else 1082 do ifft=1,nfft 1083 lambda =lambda +rhop(ifft,1)*rhor_dop_el_(ifft)*gamma(ifft,1) 1084 lambda_ipm=lambda_ipm+rhop(ifft,1)*rhor_dop_el_(ifft)*gamma(ifft,2) 1085 end do 1086 end if 1087 if (usecore==1) then 1088 do ifft=1,nfft 1089 lambda_core =lambda_core +rhop(ifft,1)*xccc3d(ifft)*gamma(ifft,1) 1090 lambda_core_ipm=lambda_core_ipm+rhop(ifft,1)*xccc3d(ifft) 1091 end do 1092 end if 1093 lambda =lambda *ucvol/dble(nfftot) 1094 lambda_ipm =lambda_ipm *ucvol/dble(nfftot) 1095 lambda_core =lambda_core *ucvol/dble(nfftot) 1096 lambda_core_ipm=lambda_core_ipm*ucvol/dble(nfftot) 1097 ABI_FREE(gamma) 1098 ABI_FREE(rhoe) 1099 ABI_FREE(rhop) 1100 if (option/=1) then 1101 ABI_FREE(rhor_dop_el_) 1102 end if 1103 ! NC pseudopotential: check electrons/positron number 1104 if (dtset%usepaw==0.and.igam==ngamma) then 1105 nbec=zero;nbev=zero;nbp=zero 1106 if (electronpositron%particle==EP_ELECTRON) then 1107 do ifft=1,nfft 1108 nbec=nbec+xccc3d(ifft) 1109 nbev=nbev+electronpositron%rhor_ep(ifft,1) 1110 nbp =nbp +rhor(ifft,1) 1111 end do 1112 else 1113 do ifft=1,nfft 1114 nbec=nbec+xccc3d(ifft) 1115 nbev=nbev+rhor(ifft,1) 1116 nbp =nbp +electronpositron%rhor_ep(ifft,1) 1117 end do 1118 end if 1119 nbec=nbec*ucvol/dble(nfftot) 1120 nbev=nbev*ucvol/dble(nfftot) 1121 nbp =nbp *ucvol/dble(nfftot) 1122 end if 1123 1124 ! MPI parallelization 1125 if(mpi_enreg%nproc_fft>1)then 1126 call xmpi_sum(lambda ,mpi_enreg%comm_fft,ierr) 1127 call xmpi_sum(lambda_ipm,mpi_enreg%comm_fft,ierr) 1128 call xmpi_sum(lambda_core ,mpi_enreg%comm_fft,ierr) 1129 call xmpi_sum(lambda_core_ipm,mpi_enreg%comm_fft,ierr) 1130 if (dtset%usepaw==0.and.igam==ngamma) then 1131 call xmpi_sum(nbec,mpi_enreg%comm_fft,ierr) 1132 call xmpi_sum(nbev,mpi_enreg%comm_fft,ierr) 1133 call xmpi_sum(nbp ,mpi_enreg%comm_fft,ierr) 1134 end if 1135 end if 1136 1137 1138 ! PAW: add on-site contributions to electron-positron annihilation rate 1139 ! ---------------------------------------------------------------------------------------- 1140 if (dtset%usepaw==1) then 1141 1142 lambda_paw =zero;lambda_paw_ipm =zero 1143 lambda_core_paw=zero;lambda_core_paw_ipm=zero 1144 1145 ! Loop on atoms 1146 do iatom=1,my_natom 1147 1148 itypat=pawrhoij(iatom)%itypat 1149 lmn2_size=pawtab(itypat)%lmn2_size 1150 mesh_size=pawtab(itypat)%mesh_size 1151 lm_size=pawtab(itypat)%lcut_size**2 1152 cplex=1 1153 ngr=0;if (ngrad==2) ngr=mesh_size 1154 1155 ! Allocations of "on-site" densities 1156 ABI_MALLOC(rho1 ,(cplex*mesh_size,lm_size,nspden_ep)) 1157 ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden_ep)) 1158 ABI_MALLOC(rho1_ep ,(cplex*mesh_size,lm_size,nspden_ep)) 1159 ABI_MALLOC(trho1_ep,(cplex*mesh_size,lm_size,nspden_ep)) 1160 if (option/=1) then 1161 ABI_MALLOC(rho1_j ,(cplex*mesh_size,lm_size,nspden_ep)) 1162 ABI_MALLOC(trho1_j,(cplex*mesh_size,lm_size,nspden_ep)) 1163 else 1164 ABI_MALLOC(rho1_j ,(0,0,0)) 1165 ABI_MALLOC(trho1_j,(0,0,0)) 1166 end if 1167 if (include_nhat_in_gamma) then 1168 ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden_ep)) 1169 ABI_MALLOC(nhat1_ep,(cplex*mesh_size,lm_size,nspden_ep)) 1170 else 1171 ABI_MALLOC(nhat1,(0,0,0)) 1172 ABI_MALLOC(nhat1_ep,(0,0,0)) 1173 end if 1174 if (include_nhat_in_gamma.and.option/=1) then 1175 ABI_MALLOC(nhat1_j,(cplex*mesh_size,lm_size,nspden_ep)) 1176 else 1177 ABI_MALLOC(nhat1_j,(0,0,0)) 1178 end if 1179 ABI_MALLOC(lmselect,(lm_size)) 1180 ABI_MALLOC(lmselect_ep,(lm_size)) 1181 ABI_MALLOC(lmselect_dum,(lm_size)) 1182 1183 ! Compute "on-site" densities (n1, ntild1, nhat1) for electron and positron ===== 1184 lmselect(:)=.true. 1185 opt_dens=1;if (include_nhat_in_gamma) opt_dens=0 1186 call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1,nspden_ep,1,& 1187 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij(iatom),& 1188 & pawtab(itypat),rho1,trho1) 1189 lmselect_ep(:)=.true. 1190 call pawdensities(rdum,cplex,iatom,lmselect_ep,lmselect_dum,lm_size,nhat1_ep,nspden_ep,1,& 1191 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_ep_(iatom),& 1192 & pawtab(itypat),rho1_ep,trho1_ep) 1193 1194 ! For state dependent scheme in Doppler ===== 1195 ! Compute "on-site" densities (n1, ntild1, nhat1) for a given electron state j===== 1196 if (option/=1) then 1197 opt_dens=1;if (include_nhat_in_gamma) opt_dens=0 1198 call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1_j,nspden_ep,1,& 1199 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_dop_el(iatom),& 1200 & pawtab(itypat),rho1_j,trho1_j) 1201 end if 1202 1203 ! Compute contribution to annihilation rate: 1204 ! Loop: first step: compute all-electron contribution (from n^1, n_c) 1205 ! 2nd step: compute pseudo contribution (from tild_n^1, hat_n^1, tild_n_c) 1206 do iloop=1,2 1207 if (iloop==1) usecore=1 1208 if (iloop==2) usecore=pawtab(itypat)%usetcore 1209 ABI_MALLOC(rhocore,(mesh_size)) 1210 1211 ! First formalism: use densities on r,theta,phi 1212 if (dtset%pawxcdev==0) then 1213 1214 ABI_MALLOC(gamma,(mesh_size,2)) 1215 ABI_MALLOC(rhoarr1,(mesh_size)) 1216 ABI_MALLOC(rhoarr1_ep,(mesh_size)) 1217 if (option/=1) then 1218 ABI_MALLOC(rhoarr1_j,(mesh_size)) 1219 end if 1220 ! Loop on the angular part 1221 do ipt=1,pawang%angl_size 1222 ! Build densities 1223 rhoarr1=zero;rhoarr1_ep=zero;rhocore=zero 1224 if (option/=1) rhoarr1_j=zero 1225 if (iloop==1) then 1226 do ilm=1,lm_size 1227 if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+rho1(:,ilm,1)*pawang%ylmr(ilm,ipt) 1228 end do 1229 if (option/=1) then 1230 do ilm=1,lm_size 1231 if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1_j(:)+rho1_j(:,ilm,1)*pawang%ylmr(ilm,ipt) 1232 end do 1233 end if 1234 do ilm=1,lm_size 1235 if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+rho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt) 1236 end do 1237 if (usecore==1) rhocore(:)=pawtab(itypat)%coredens(:) 1238 else 1239 if (include_nhat_in_gamma) then 1240 do ilm=1,lm_size 1241 if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+(trho1(:,ilm,1)+nhat1(:,ilm,1))*pawang%ylmr(ilm,ipt) 1242 end do 1243 if (option/=1) then 1244 do ilm=1,lm_size 1245 if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1_j(:)+(trho1_j(:,ilm,1)+nhat1_j(:,ilm,1))*pawang%ylmr(ilm,ipt) 1246 end do 1247 end if 1248 do ilm=1,lm_size 1249 if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+(trho1_ep(:,ilm,1)+nhat1_ep(:,ilm,1))*pawang%ylmr(ilm,ipt) 1250 end do 1251 else 1252 do ilm=1,lm_size 1253 if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+trho1(:,ilm,1)*pawang%ylmr(ilm,ipt) 1254 end do 1255 if (option/=1) then 1256 do ilm=1,lm_size 1257 if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1(:)+trho1_j(:,ilm,1)*pawang%ylmr(ilm,ipt) 1258 end do 1259 end if 1260 do ilm=1,lm_size 1261 if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+trho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt) 1262 end do 1263 end if 1264 if (usecore==1) rhocore(:)=pawtab(itypat)%tcoredens(:,1) 1265 end if 1266 ! Make the densities positive 1267 if (electronpositron%particle==EP_ELECTRON) then 1268 call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos) 1269 call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos) 1270 else if (electronpositron%particle==EP_POSITRON) then 1271 call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos) 1272 call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos) 1273 if (option/=1) then 1274 call mkdenpos(iwarnj,mesh_size,1,1,rhoarr1_j,dtset%xc_denpos) 1275 end if 1276 end if 1277 ! Compute Gamma 1278 ABI_MALLOC(grhoe2,(ngr)) 1279 ABI_MALLOC(grhocore2,(ngr)) 1280 if (option==1.or.option==2) then 1281 if (electronpositron%particle==EP_ELECTRON) then 1282 call gammapositron(gamma,grhocore2,grhoe2,igamma(igam),ngr,mesh_size,& 1283 & rhocore,rhoarr1_ep,rhoarr1,usecore) 1284 else if (electronpositron%particle==EP_POSITRON) then 1285 call gammapositron(gamma,grhocore2,grhoe2,igamma(igam),ngr,mesh_size,& 1286 & rhocore,rhoarr1,rhoarr1_ep,usecore) 1287 end if 1288 else 1289 gamma(:,:)=one 1290 end if 1291 ABI_FREE(grhoe2) 1292 ABI_FREE(grhocore2) 1293 ! Compute contribution to annihilation rates 1294 ABI_MALLOC(ff,(mesh_size)) 1295 if (option/=1) rhoarr1(:)=rhoarr1_j(:) 1296 do ii=1,4 1297 if (ii==1) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhoarr1_ep(1:mesh_size) & 1298 & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 1299 if (ii==2) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhoarr1_ep(1:mesh_size) & 1300 & *gamma(1:mesh_size,2)*pawrad(itypat)%rad(1:mesh_size)**2 1301 if (electronpositron%particle==EP_ELECTRON) then 1302 if (ii==3) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhocore(1:mesh_size) & 1303 & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 1304 if (ii==4) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhocore(1:mesh_size) & 1305 & *pawrad(itypat)%rad(1:mesh_size)**2 1306 else 1307 if (ii==3) ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocore(1:mesh_size) & 1308 & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 1309 if (ii==4) ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocore(1:mesh_size) & 1310 & *pawrad(itypat)%rad(1:mesh_size)**2 1311 end if 1312 call simp_gen(intg,ff,pawrad(itypat)) 1313 intg=intg*pawang%angwgth(ipt)*four_pi 1314 if (ii==1) lambda_paw =lambda_paw +lsign(iloop)*intg 1315 if (ii==2) lambda_paw_ipm =lambda_paw_ipm +lsign(iloop)*intg 1316 if (ii==3) lambda_core_paw =lambda_core_paw +lsign(iloop)*intg 1317 if (ii==4) lambda_core_paw_ipm=lambda_core_paw_ipm+lsign(iloop)*intg 1318 end do 1319 ABI_FREE(ff) 1320 end do ! ipt 1321 ABI_FREE(gamma) 1322 ABI_FREE(rhoarr1) 1323 ABI_FREE(rhoarr1_ep) 1324 if (option/=1) then 1325 ABI_FREE(rhoarr1_j) 1326 end if 1327 1328 ! Second formalism: use (l,m) moments for densities 1329 else if (dtset%pawxcdev/=0) then 1330 1331 ! Build densities 1332 ABI_MALLOC(gammam,(mesh_size,2,lm_size)) 1333 ABI_MALLOC(rhotot,(mesh_size,lm_size)) 1334 ABI_MALLOC(rhotot_ep,(mesh_size,lm_size)) 1335 ABI_MALLOC(rhosph,(mesh_size)) 1336 ABI_MALLOC(rhosph_ep,(mesh_size)) 1337 if (option/=1) then 1338 ABI_MALLOC(rhotot_j,(mesh_size,lm_size)) 1339 ABI_MALLOC(rhosph_j,(mesh_size)) 1340 end if 1341 if (usecore==0) rhocore(:)=zero 1342 if (iloop==1) then 1343 rhotot (:,:)=rho1 (:,:,1) 1344 rhotot_ep(:,:)=rho1_ep(:,:,1) 1345 if (option/=1) rhotot_j (:,:)=rho1_j (:,:,1) 1346 if (usecore==1) rhocore(:)=pawtab(itypat)%coredens(:) 1347 else 1348 if (include_nhat_in_gamma) then 1349 rhotot (:,:)=trho1 (:,:,1)+nhat1 (:,:,1) 1350 rhotot_ep(:,:)=trho1_ep(:,:,1)+nhat1_ep(:,:,1) 1351 if (option/=1) rhotot_j (:,:)=trho1_j (:,:,1)+nhat1_j (:,:,1) 1352 else 1353 rhotot (:,:)=trho1 (:,:,1) 1354 rhotot_ep(:,:)=trho1_ep(:,:,1) 1355 if (option/=1) rhotot_j (:,:)=trho1_j (:,:,1) 1356 end if 1357 if (usecore==1) rhocore(:)=pawtab(itypat)%tcoredens(:,1) 1358 end if 1359 rhosph (:)=rhotot (:,1)/sqfpi 1360 rhosph_ep(:)=rhotot_ep(:,1)/sqfpi 1361 if (option/=1) rhosph_j (:)=rhotot_j(:,1)/sqfpi 1362 ! Make spherical densities positive 1363 if (electronpositron%particle==EP_ELECTRON) then 1364 call mkdenpos(iwarnp,mesh_size,1,1,rhosph ,dtset%xc_denpos) 1365 call mkdenpos(iwarn ,mesh_size,1,1,rhosph_ep,dtset%xc_denpos) 1366 else if (electronpositron%particle==EP_POSITRON) then 1367 call mkdenpos(iwarn ,mesh_size,1,1,rhosph ,dtset%xc_denpos) 1368 call mkdenpos(iwarnp,mesh_size,1,1,rhosph_ep,dtset%xc_denpos) 1369 if (option/=1) then 1370 call mkdenpos(iwarnp,mesh_size,1,1,rhosph_j,dtset%xc_denpos) 1371 end if 1372 end if 1373 ! Need gradients of electronic densities for GGA 1374 ABI_MALLOC(grhoe2,(ngr)) 1375 ABI_MALLOC(grhocore2,(ngr)) 1376 if (ngr>0) then 1377 if (electronpositron%particle==EP_ELECTRON) then 1378 call nderiv_gen(grhoe2,rhosph_ep,pawrad(itypat)) 1379 else if (electronpositron%particle==EP_POSITRON) then 1380 call nderiv_gen(grhoe2,rhosph,pawrad(itypat)) 1381 end if 1382 grhoe2(:)=grhoe2(:)**2 1383 if (usecore==1) then 1384 call nderiv_gen(grhocore2,rhocore,pawrad(itypat)) 1385 grhocore2(:)=grhocore2(:)**2 1386 end if 1387 end if 1388 ! Compute Gamma for (rho-,rho+), 1389 ! (rho- +drho-,rho+), (rho- -drho-,rho+), 1390 ! (rho-,rho+ +drho+), (rho-,rho+ -drho+), 1391 ! (rho- +drho-,rho+ +drho+), (rho- -drho-,rho+ -drho+) 1392 ! Do a seven steps loop 1393 ABI_MALLOC(gam_,(mesh_size,2,7)) 1394 ABI_MALLOC(rho_,(mesh_size)) 1395 ABI_MALLOC(rho_ep_,(mesh_size)) 1396 ABI_MALLOC(rhocor_,(mesh_size)) 1397 ABI_MALLOC(grho2_,(ngr)) 1398 ABI_MALLOC(grhocor2_,(ngr)) 1399 do ii=1,7 1400 ! Apply delta to get perturbed densities 1401 rho_(:)=rhosph(:);rho_ep_(:)=rhosph_ep(:);if (usecore==1) rhocor_(:)=rhocore(:) 1402 if (ngr>0) grho2_(:)=grhoe2(:) 1403 if (ngr>0) grhocor2_(:)=grhocore2(:) 1404 if (ii==2.or.ii==4.or.ii==6) fact=(one+delta) 1405 if (ii==3.or.ii==5.or.ii==7) fact=(one-delta) 1406 fact2=fact**2 1407 if (ii==2.or.ii==3.or.ii==6.or.ii==7) then 1408 rho_(:)=fact*rho_(:) 1409 if (electronpositron%particle==EP_POSITRON) then 1410 if (ngr>0) grho2_(:)=fact2*grho2_(:) 1411 if (usecore==1)rhocor_(:)=fact*rhocor_(:) 1412 if (ngr>0.and.usecore==1) grhocor2_(:)=fact2*grhocor2_(:) 1413 end if 1414 end if 1415 if (ii==4.or.ii==5.or.ii==6.or.ii==7) then 1416 rho_ep_(:)=fact*rho_ep_(:) 1417 if (electronpositron%particle==EP_ELECTRON) then 1418 if (ngr>0) grho2_(:)=fact2*grho2_(:) 1419 if (usecore==1)rhocor_(:)=fact*rhocor_(:) 1420 if (ngr>0.and.usecore==1) grhocor2_(:)=fact2*grhocor2_(:) 1421 end if 1422 end if 1423 ! Compute gamma for these perturbed densities 1424 if (option==1.or.option==2) then 1425 if (electronpositron%particle==EP_ELECTRON) then 1426 call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma(igam),ngr,mesh_size,rhocor_,rho_ep_,rho_,usecore) 1427 else if (electronpositron%particle==EP_POSITRON) then 1428 call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma(igam),ngr,mesh_size,rhocor_,rho_,rho_ep_,usecore) 1429 end if 1430 else 1431 gam_(:,:,:)=one 1432 end if 1433 end do ! end loop ii=1,7 1434 1435 ABI_FREE(rhocor_) 1436 ABI_FREE(grho2_) 1437 ABI_FREE(grhocor2_) 1438 ABI_FREE(grhoe2) 1439 ABI_FREE(grhocore2) 1440 rho_ (:)=rhosph (:);if (electronpositron%particle==EP_POSITRON.and.usecore==1) rho_ (:)=rho_ (:)+rhocore(:) 1441 rho_ep_(:)=rhosph_ep(:);if (electronpositron%particle==EP_ELECTRON.and.usecore==1) rho_ep_(:)=rho_ep_(:)+rhocore(:) 1442 ! Compute numerical first and second derivatives of Gamma 1443 ! d1gam(1) = dgam/drho+ (particle=ELECTRON), dgam/drho- (particle=POSITRON) 1444 ! d1gam(2) = dgam/drho- (particle=ELECTRON), dgam/drho+ (particle=POSITRON) 1445 ABI_MALLOC(d1gam,(mesh_size,2,2)) 1446 d1gam(:,:,:)=zero 1447 do ir=1,mesh_size 1448 if (rho_ (ir)>tol14) d1gam(ir,1,1)=(gam_(ir,1,2)-gam_(ir,1,3))*half/(delta*rho_ (ir)) 1449 if (rhosph (ir)>tol14) d1gam(ir,2,1)=(gam_(ir,2,2)-gam_(ir,2,3))*half/(delta*rhosph (ir)) 1450 if (rho_ep_ (ir)>tol14) d1gam(ir,1,2)=(gam_(ir,1,4)-gam_(ir,1,5))*half/(delta*rho_ep_ (ir)) 1451 if (rhosph_ep(ir)>tol14) d1gam(ir,2,2)=(gam_(ir,2,4)-gam_(ir,2,5))*half/(delta*rhosph_ep(ir)) 1452 end do 1453 1454 ! d2gam(1) = d2gam/drho+_drho+ (particle=ELECTRON), dgam/drho-_drho- (particle=POSITRON) 1455 ! d2gam(2) = d2gam/drho-_drho+ (particle=ELECTRON), dgam/drho+_drho- (particle=POSITRON) 1456 ! d2gam(3) = d2gam/drho-_drho- (particle=ELECTRON), dgam/drho+_drho+ (particle=POSITRON) 1457 ABI_MALLOC(d2gam,(mesh_size,2,3)) 1458 d2gam(:,:,:)=zero 1459 do ir=1,mesh_size 1460 if (rho_ (ir)>tol14) d2gam(ir,1,1)=(gam_(ir,1,2)+gam_(ir,1,3)-two*gam_(ir,1,1))/(delta*rho_ (ir))**2 1461 if (rhosph(ir)>tol14) d2gam(ir,2,1)=(gam_(ir,2,2)+gam_(ir,2,3)-two*gam_(ir,2,1))/(delta*rhosph(ir))**2 1462 if (rho_ep_(ir)>tol14) then 1463 d2gam(ir,1,3)=(gam_(ir,1,4)+gam_(ir,1,5)-two*gam_(ir,1,1))/(delta*rho_ep_(ir))**2 1464 if (rho_(ir)>tol14) then 1465 d2gam(ir,1,2)=(gam_(ir,1,6)+gam_(ir,1,7)+two*gam_(ir,1,1) & 1466 & -gam_(ir,1,2)-gam_(ir,1,3)-gam_(ir,1,4)-gam_(ir,1,5)) & 1467 & *half/(delta*rho_(ir))/(delta*rho_ep_(ir)) 1468 end if 1469 end if 1470 if (rhosph_ep(ir)>tol14) then 1471 d2gam(ir,2,3)=(gam_(ir,2,4)+gam_(ir,2,5)-two*gam_(ir,2,1))/(delta*rhosph_ep(ir))**2 1472 if (rhosph(ir)>tol14) then 1473 d2gam(ir,2,2)=(gam_(ir,2,6)+gam_(ir,2,7)+two*gam_(ir,2,1) & 1474 & -gam_(ir,2,2)-gam_(ir,2,3)-gam_(ir,2,4)-gam_(ir,2,5)) & 1475 & *half/(delta*rhosph(ir))/(delta*rhosph_ep(ir)) 1476 end if 1477 end if 1478 end do 1479 ABI_FREE(rho_) 1480 ABI_FREE(rho_ep_) 1481 ! Compute useful sums of densities 1482 ABI_MALLOC(v1sum,(mesh_size,3)) 1483 if ( dtset%pawxcdev>=2) then 1484 ABI_MALLOC(v2sum,(mesh_size,lm_size,3)) 1485 else 1486 ABI_MALLOC(v2sum,(0,0,0)) 1487 end if 1488 rhotot(:,1)=sqfpi*rhosph(:);rhotot_ep(:,1)=sqfpi*rhosph_ep(:) 1489 call pawxcsum(1,1,1,lmselect,lmselect_ep,lm_size,mesh_size,3,dtset%pawxcdev,& 1490 & pawang,rhotot,rhotot_ep,v1sum,v2sum) 1491 ! Compute final developpment of gamma moments 1492 gammam(:,:,:)=zero 1493 gammam(:,:,1)=gam_(:,:,1)*sqfpi 1494 gammam(:,1,1)=gammam(:,1,1)+(d2gam(:,1,2)*v1sum(:,2) & 1495 & +half*(d2gam(:,1,1)*v1sum(:,1)+d2gam(:,1,3)*v1sum(:,3)))/sqfpi 1496 gammam(:,2,1)=gammam(:,2,1)+(d2gam(:,2,2)*v1sum(:,2) & 1497 & +half*(d2gam(:,2,1)*v1sum(:,1)+d2gam(:,2,3)*v1sum(:,3)))/sqfpi 1498 do ilm=2,lm_size 1499 if (lmselect(ilm)) then 1500 gammam(:,1,ilm)=gammam(:,1,ilm)+d1gam(:,1,1)*rhotot(:,ilm) 1501 gammam(:,2,ilm)=gammam(:,2,ilm)+d1gam(:,2,1)*rhotot(:,ilm) 1502 end if 1503 if (lmselect_ep(ilm)) then 1504 gammam(:,1,ilm)=gammam(:,1,ilm)+d1gam(:,1,2)*rhotot_ep(:,ilm) 1505 gammam(:,2,ilm)=gammam(:,2,ilm)+d1gam(:,2,2)*rhotot_ep(:,ilm) 1506 end if 1507 end do 1508 if (dtset%pawxcdev>1) then 1509 do ilm=2,lm_size 1510 gammam(:,1,ilm)=gammam(:,1,ilm)+d2gam(:,1,2)*v2sum(:,ilm,2) & 1511 & +half*(d2gam(:,1,1)*v2sum(:,ilm,1)+d2gam(:,1,3)*v2sum(:,ilm,3)) 1512 gammam(:,2,ilm)=gammam(:,2,ilm)+d2gam(:,2,2)*v2sum(:,ilm,2) & 1513 & +half*(d2gam(:,2,1)*v2sum(:,ilm,1)+d2gam(:,2,3)*v2sum(:,ilm,3)) 1514 end do 1515 end if 1516 ABI_FREE(gam_) 1517 ABI_FREE(d1gam) 1518 ABI_FREE(d2gam) 1519 ABI_FREE(v1sum) 1520 ABI_FREE(v2sum) 1521 1522 ! Compute contribution to annihilation rate 1523 1524 ! In state dependent scheme for Doppler, replace electronic density with one state 1525 if (option/=1) then 1526 rhotot(:,:) = rhotot_j(:,:) 1527 rhosph (:) = rhosph_j (:) 1528 end if 1529 1530 ABI_MALLOC(gg,(mesh_size,4)) 1531 gg=zero 1532 ABI_MALLOC(rhoarr1,(mesh_size)) 1533 ABI_MALLOC(rhoarr2,(mesh_size)) 1534 do ilm=1,lm_size 1535 do ilm1=1,lm_size 1536 if (lmselect(ilm1)) then 1537 if (ilm1==1) rhoarr1(:)=sqfpi*rhosph(:) 1538 if (ilm1/=1) rhoarr1(:)=rhotot(:,ilm1) 1539 do ilm2=1,lm_size 1540 if (lmselect_ep(ilm2)) then 1541 if (ilm2==1) rhoarr2(:)=sqfpi*rhosph_ep(:) 1542 if (ilm2/=1) rhoarr2(:)=rhotot_ep(:,ilm2) 1543 if (ilm1>=ilm2) then 1544 isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2) 1545 else 1546 isel=pawang%gntselect(ilm,ilm1+ilm2*(ilm2-1)/2) 1547 end if 1548 if (isel>0) then 1549 fact=pawang%realgnt(isel) 1550 gg(:,1)=gg(:,1)+fact*rhoarr1(:)*rhoarr2(:)*gammam(:,1,ilm) 1551 gg(:,2)=gg(:,2)+fact*rhoarr1(:)*rhoarr2(:)*gammam(:,2,ilm) 1552 end if 1553 end if 1554 end do 1555 end if 1556 end do 1557 end do 1558 ABI_FREE(rhoarr1) 1559 ABI_FREE(rhoarr2) 1560 if (electronpositron%particle==EP_ELECTRON) then 1561 do ilm=1,lm_size 1562 if (lmselect(ilm)) gg(:,3)=gg(:,3)+rhotot(:,ilm)*rhocore(:)*gammam(:,1,ilm) 1563 end do 1564 gg(:,4)=sqfpi*rhotot(:,1)*rhocore(:) 1565 else if (electronpositron%particle==EP_POSITRON) then 1566 do ilm=1,lm_size 1567 if (lmselect_ep(ilm)) gg(:,3)=gg(:,3)+rhotot_ep(:,ilm)*rhocore(:)*gammam(:,1,ilm) 1568 end do 1569 gg(:,4)=sqfpi*rhotot_ep(:,1)*rhocore(:) 1570 end if 1571 do ii=1,4 1572 gg(1:mesh_size,ii)=gg(1:mesh_size,ii)*pawrad(itypat)%rad(1:mesh_size)**2 1573 call simp_gen(intg,gg(:,ii),pawrad(itypat)) 1574 if (ii==1) lambda_paw =lambda_paw +lsign(iloop)*intg 1575 if (ii==2) lambda_paw_ipm =lambda_paw_ipm +lsign(iloop)*intg 1576 if (ii==3) lambda_core_paw =lambda_core_paw +lsign(iloop)*intg 1577 if (ii==4) lambda_core_paw_ipm=lambda_core_paw_ipm+lsign(iloop)*intg 1578 end do 1579 ABI_FREE(gg) 1580 ABI_FREE(gammam) 1581 ABI_FREE(rhotot) 1582 ABI_FREE(rhotot_ep) 1583 ABI_FREE(rhosph) 1584 ABI_FREE(rhosph_ep) 1585 if (option/=1) then 1586 ABI_FREE(rhotot_j) 1587 ABI_FREE(rhosph_j) 1588 end if 1589 1590 end if ! dtset%pawxcdev 1591 1592 ABI_FREE(rhocore) 1593 1594 end do ! iloop 1595 1596 ABI_FREE(rho1) 1597 ABI_FREE(trho1) 1598 ABI_FREE(rho1_ep) 1599 ABI_FREE(trho1_ep) 1600 ABI_FREE(rho1_j) 1601 ABI_FREE(trho1_j) 1602 ABI_FREE(nhat1) 1603 ABI_FREE(nhat1_ep) 1604 ABI_FREE(nhat1_j) 1605 ABI_FREE(lmselect) 1606 ABI_FREE(lmselect_ep) 1607 ABI_FREE(lmselect_dum) 1608 1609 end do ! iatom 1610 1611 ! Reduction in case of distribution over atomic sites 1612 if (mpi_enreg%nproc_atom>1) then 1613 mpibuf(1)=lambda_paw ;mpibuf(2)=lambda_paw_ipm 1614 mpibuf(3)=lambda_core_paw;mpibuf(4)=lambda_core_paw_ipm 1615 call xmpi_sum(mpibuf,mpi_enreg%comm_atom,ierr) 1616 lambda_paw=mpibuf(1) ;lambda_paw_ipm=mpibuf(2) 1617 lambda_core_paw=mpibuf(3);lambda_core_paw_ipm=mpibuf(4) 1618 end if 1619 1620 ! Add plane-wave and PAW contributions to annihilation rates 1621 lambda =lambda +lambda_paw 1622 lambda_ipm =lambda_ipm +lambda_paw_ipm 1623 lambda_core =lambda_core +lambda_core_paw 1624 lambda_core_ipm=lambda_core_ipm+lambda_core_paw_ipm 1625 end if ! dtset%usepaw 1626 1627 1628 ! Convert into proper units and print 1629 ! --------------------------------------------------------------------------------------- 1630 1631 ! Sum valence and core contributions to annihilation rates 1632 lambda =lambda +lambda_core 1633 lambda_ipm =lambda_ipm +lambda_core_ipm 1634 if (dtset%usepaw==1) then 1635 lambda_paw =lambda_paw +lambda_core_paw 1636 lambda_paw_ipm=lambda_paw_ipm+lambda_core_paw_ipm 1637 end if 1638 1639 ! Set annihilation rate in proper unit (picosec.) 1640 units=pi*(one/InvFineStruct)**3/Time_Sec/1.e12_dp/electronpositron%posocc 1641 1642 lambda =lambda *units 1643 lambda_ipm =lambda_ipm *units 1644 lambda_core =lambda_core *units 1645 lambda_core_ipm=lambda_core_ipm*units 1646 lifetime =one/lambda 1647 lifetime_ipm=one/lambda_ipm 1648 electronpositron%lambda=lambda 1649 electronpositron%lifetime=lifetime 1650 if (dtset%usepaw==1) then 1651 lambda_paw =lambda_paw *units 1652 lambda_paw_ipm =lambda_paw_ipm *units 1653 lambda_core_paw =lambda_core_paw *units 1654 lambda_core_paw_ipm=lambda_core_paw_ipm*units 1655 end if 1656 rate=lambda-lambda_core-lambda_paw+lambda_core_paw 1657 rate_paw=lambda_paw-lambda_core_paw 1658 ! Print life time and additional information 1659 if (option==1) then 1660 if (igam==1) then 1661 write(msg,'(a,80("-"),2a)') ch10,ch10,' Results for electron-positron annihilation:' 1662 call wrtout(ab_out,msg,'COLL') 1663 call wrtout(std_out,msg,'COLL') 1664 end if 1665 if (ngamma>1.and.igam==1) then 1666 write(msg,'(a,i1,a)') ch10,ngamma,& 1667 & ' computations of positron lifetime have been performed (with different enhancement factors).' 1668 call wrtout(ab_out,msg,'COLL') 1669 call wrtout(std_out,msg,'COLL') 1670 end if 1671 if (ngamma>1) then 1672 write(msg,'(2a,i1)') ch10,"########## Lifetime computation ",igam 1673 call wrtout(ab_out,msg,'COLL') 1674 call wrtout(std_out,msg,'COLL') 1675 end if 1676 if (abs(electronpositron%ixcpositron)==1) then 1677 write(msg,'(4a)') ch10,' # Zero-positron density limit of Arponen and Pajanne provided by Boronski & Nieminen',& 1678 & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]] 1679 else if (electronpositron%ixcpositron==11) then 1680 write(msg,'(4a)') ch10,' # Zero-positron density limit of Arponen and Pajanne fitted by Sterne & Kaiser',& 1681 & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)' ! [[cite:Sterne1991]] 1682 else if (electronpositron%ixcpositron==2) then 1683 write(msg,'(4a)') ch10,' # Electron-positron correlation provided by Puska, Seitsonen, and Nieminen',& 1684 & ch10,' Ref: M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994)' ! [[cite:Puska1994]] 1685 else if (electronpositron%ixcpositron==3) then 1686 write(msg,'(8a)') ch10,' # Zero-positron density limit of Arponen and Pajanne provided by Boronski & Nieminen',& 1687 & ch10,' + GGA corrections',& 1688 & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)',& ! [[cite:Boronski1986]] 1689 & ch10,' B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1995)' ! [[cite:Barbiellini1995]] 1690 else if (electronpositron%ixcpositron==31) then 1691 write(msg,'(8a)') ch10,' # Zero-positron density limit of Arponen and Pajanne fitted by Sterne & Kaiser',& 1692 & ch10,' + GGA corrections',& 1693 & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)',& ! [[cite:Sterne1991]] 1694 & ch10,' B. Barbiellini, M.J. Puska, T. Torsti and R.M. Nieminen, Phys. Rev. B 51, 7341 (1995)' ! [[cite:Barbiellini1995]] 1695 end if 1696 call wrtout(ab_out,msg,'COLL') 1697 call wrtout(std_out, msg,'COLL') 1698 if (igamma(igam)==0) then 1699 write(msg,'(a)') ' # Enhancement factor set to one (test)' 1700 else if (igamma(igam)==1) then 1701 write(msg,'(3a)') ' # Enhancement factor of Boronski & Nieminen',& 1702 & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]] 1703 else if (igamma(igam)==2) then 1704 write(msg,'(3a)') ' # Enhancement factor of Boronski & Nieminen IN THE RPA LIMIT',& 1705 & ch10,' Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]] 1706 else if (igamma(igam)==3) then 1707 write(msg,'(3a)') ' # Enhancement factor of Sterne & Kaiser',& 1708 & ch10,' Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)' ! [[cite:Sterne1991]] 1709 else if (igamma(igam)==4) then 1710 write(msg,'(3a)') ' # Enhancement factor of Puska, Seitsonen, and Nieminen',& 1711 & ch10,' Ref.: M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994)' ! [[cite:Puska1994]] 1712 end if 1713 call wrtout(ab_out,msg,'COLL') 1714 call wrtout(std_out,msg,'COLL') 1715 write(msg, '(4(2a,es16.8))' ) ch10,& 1716 & ' Positron lifetime (ps) =',lifetime ,ch10,& 1717 & ' Positron lifetime with IPM for core elec. (ps) =',lifetime_ipm,ch10,& 1718 & ' Annihilation rate (ns-1) =',lambda *1000._dp,ch10,& 1719 & ' Annihilation rate with IPM for core elec. (ns-1) =',lambda_ipm*1000._dp 1720 call wrtout(ab_out,msg,'COLL') 1721 call wrtout(std_out,msg,'COLL') 1722 write(msg,'(2a,5(2a,es16.8))' ) ch10,& 1723 & ' Annihilation rate core/valence decomposition:',ch10,& 1724 & ' Core contribution to ann.rate (ns-1) =', lambda_core *1000._dp,ch10,& 1725 & ' Valence contribution to ann.rate (ns-1) =',(lambda-lambda_core) *1000._dp,ch10,& 1726 & ' Core contribution to ann.rate with IPM (ns-1) =', lambda_core_ipm *1000._dp,ch10,& 1727 & ' Valence contribution to ann.rate with IPM (ns-1) =',(lambda_ipm-lambda_core_ipm) *1000._dp 1728 call wrtout(ab_out,msg,'COLL') 1729 call wrtout(std_out,msg,'COLL') 1730 if (dtset%usepaw==1) then 1731 write(msg, '(2a,6(2a,es16.8))' ) ch10,& 1732 & ' Annihilation rate PAW decomposition:',ch10,& 1733 & ' Plane-wave contribution to ann.rate (ns-1) =',(lambda-lambda_paw)*1000._dp,ch10,& 1734 & ' Plane-wave valence contribution to ann.rate (ns-1) =',(lambda-lambda_paw-lambda_core+lambda_core_paw)*1000._dp,ch10,& 1735 & ' On-site core contribution to ann.rate (ns-1) =', lambda_core_paw*1000._dp,ch10,& 1736 & ' On-site valence contribution to ann.rate (ns-1) =',(lambda_paw-lambda_core_paw)*1000._dp,ch10,& 1737 & ' Plane-wave contribution to ann.rate with IPM (ns-1) =',(lambda_ipm-lambda_paw_ipm)*1000._dp,ch10,& 1738 & ' Plane-wave core contrb. to ann.rate with IPM (ns-1) =',(lambda_core_ipm-lambda_core_paw_ipm)*1000._dp 1739 call wrtout(ab_out,msg,'COLL') 1740 call wrtout(std_out,msg,'COLL') 1741 end if 1742 if (dtset%usepaw==0.and.igam==ngamma) then ! These tests are not relevant with PAW 1743 write(msg, '(2a,3(2a,es16.8))' ) ch10,& 1744 & ' ########## Some checks, for testing purpose:',ch10,& 1745 & ' Number of core electrons =',nbec,ch10,& 1746 & ' Number of valence electrons =',nbev,ch10,& 1747 & ' Number of positrons =',nbp 1748 call wrtout(ab_out,msg,'COLL') 1749 call wrtout(std_out,msg,'COLL') 1750 end if 1751 end if !end if option 1752 end do ! Big loop on igam 1753 1754 if (option==1) then 1755 write(msg, '(3a)' ) ch10,' (*) IPM=Independent particle Model',ch10 1756 call wrtout(ab_out,msg,'COLL') 1757 call wrtout(std_out,msg,'COLL') 1758 end if !end if option 1759 1760 !Deallocate memory 1761 ABI_FREE(igamma) 1762 if (dtset%usepaw==1.and.(.not.include_nhat_in_gamma)) then 1763 ABI_FREE(rhor_) 1764 ABI_FREE(rhor_ep_) 1765 end if 1766 1767 DBG_EXIT("COLL") 1768 1769 end subroutine poslifetime
ABINIT/posratecore [ Functions ]
NAME
posratecore
FUNCTION
Calculate the annihilataion rate of a given core state
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset | nspden=number of spin-density components | ntypat=number of atom types | paral_kgb=flag controlling (k,g,bands) parallelization | pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) | usepaw=flag for PAW iatom= index of the current atom in posdoppler mesh_sizej= size of the radial mesh for the current atom in posdoppler mpi_enreg= information about MPI parallelization my_natom=number of atoms treated by current processor option= if 1, use gamma if 2, use IPM (gamma=1) pawang <type(pawang)>=paw angular mesh and related data pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
rate= annihilation rate of a given core state needed for state dependent scheme for doppler broadening
SIDE EFFECTS
SOURCE
3333 subroutine posratecore(dtset,electronpositron,iatom,my_natom,mesh_sizej,mpi_enreg,& 3334 & option,pawang,pawrad,pawrhoij,pawrhoij_ep,& 3335 & pawtab,rate,rhocorej) 3336 3337 !Arguments ------------------------------------ 3338 !scalars 3339 integer,intent(in) :: iatom,my_natom,option,mesh_sizej 3340 real(dp),intent(out) :: rate 3341 type(dataset_type), intent(in) :: dtset 3342 type(electronpositron_type),pointer :: electronpositron 3343 type(MPI_type),intent(in) :: mpi_enreg 3344 type(pawang_type), intent(in) :: pawang 3345 !arrays 3346 real(dp),intent(in) :: rhocorej(mesh_sizej) 3347 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw) 3348 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*dtset%usepaw) 3349 type(pawrhoij_type),intent(in),target :: pawrhoij_ep(my_natom*dtset%usepaw) 3350 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 3351 3352 !Local variables------------------------------- 3353 !scalars 3354 integer :: cplex,ierr,igamma,ii,ilm,ipt,ir 3355 integer :: itypat,iwarn,iwarnj,iwarnp,lm_size,lmn2_size,mesh_size 3356 integer :: ngr,ngrad,nspden_ep,opt_dens 3357 logical,parameter :: include_nhat_in_gamma=.false. 3358 real(dp),parameter :: delta=1.d-4 3359 real(dp) :: fact,fact2,intg 3360 real(dp) :: mpibuf,rdum,sqfpi 3361 character(len=500) :: msg 3362 !arrays 3363 logical,allocatable :: lmselect(:),lmselect_ep(:),lmselect_dum(:) 3364 real(dp),parameter :: qphon(3)=(/zero,zero,zero/),lsign(2)=(/one,-one/) 3365 real(dp),allocatable :: d1gam(:,:),d2gam(:,:),ff(:),gam_(:,:,:),gamma(:,:),gammam(:,:),gg(:,:) 3366 real(dp),allocatable :: grhocore2(:),grhocor2_(:),grhoe2(:),grho2_(:) 3367 real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:) 3368 real(dp),allocatable :: rho_(:),rho_ep_(:),rho1(:,:,:),rho1_ep(:,:,:) 3369 real(dp),allocatable :: rhoarr1(:),rhoarr1_ep(:),rhoarr2(:) 3370 real(dp),allocatable :: rhocore(:),rhocor_(:) 3371 real(dp),allocatable :: rhosph(:),rhosph_ep(:),rhotot(:,:),rhotot_ep(:,:) 3372 real(dp),allocatable :: trho1(:,:,:),trho1_ep(:,:,:) 3373 real(dp),allocatable :: v1sum(:,:),v2sum(:,:,:) 3374 type(pawrhoij_type),pointer :: pawrhoij_ep_(:) 3375 3376 ! ************************************************************************* 3377 3378 DBG_ENTER("COLL") 3379 3380 !Tests for developers 3381 if (.not.associated(electronpositron)) then 3382 msg='electronpositron variable must be associated!' 3383 ABI_BUG(msg) 3384 end if 3385 !Constants 3386 fact=0.0 3387 cplex=1;nspden_ep=1 3388 ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2 3389 iwarn=0;iwarnj=0;iwarnp=1 3390 sqfpi=sqrt(four_pi) 3391 3392 !Compatibility tests 3393 if (electronpositron%particle==EP_NOTHING) then 3394 msg='Not valid for electronpositron%particle=NOTHING!' 3395 ABI_BUG(msg) 3396 end if 3397 3398 if (dtset%usepaw==1) then 3399 if(dtset%pawxcdev==0.and.ngrad==2) then 3400 msg='GGA is not implemented for pawxcdev=0 (use dtset%pawxcdev/=0)!' 3401 ABI_BUG(msg) 3402 end if 3403 end if 3404 3405 !Select type(s) of enhancement factor 3406 if (electronpositron%ixcpositron==-1) igamma=0 3407 if (electronpositron%ixcpositron== 2) igamma=4 3408 if (electronpositron%ixcpositron==11.or.electronpositron%ixcpositron==31) igamma=3 3409 if (electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3) igamma=2 3410 if (option==2) igamma=0 3411 3412 pawrhoij_ep_ => pawrhoij_ep 3413 3414 rate=zero 3415 3416 itypat=pawrhoij(iatom)%itypat 3417 lmn2_size=pawtab(itypat)%lmn2_size 3418 mesh_size=pawtab(itypat)%mesh_size 3419 lm_size=pawtab(itypat)%lcut_size**2 3420 cplex=1 3421 ngr=0;if (ngrad==2) ngr=mesh_size 3422 3423 !Allocations of "on-site" densities 3424 ABI_MALLOC(rho1 ,(cplex*mesh_size,lm_size,nspden_ep)) 3425 ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden_ep)) 3426 ABI_MALLOC(rho1_ep ,(cplex*mesh_size,lm_size,nspden_ep)) 3427 ABI_MALLOC(trho1_ep,(cplex*mesh_size,lm_size,nspden_ep)) 3428 ABI_MALLOC(lmselect,(lm_size)) 3429 ABI_MALLOC(lmselect_ep,(lm_size)) 3430 ABI_MALLOC(lmselect_dum,(lm_size)) 3431 if (include_nhat_in_gamma) then 3432 ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden_ep)) 3433 ABI_MALLOC(nhat1_ep,(cplex*mesh_size,lm_size,nspden_ep)) 3434 else 3435 ABI_MALLOC(nhat1,(0,0,0)) 3436 ABI_MALLOC(nhat1_ep,(0,0,0)) 3437 end if 3438 3439 !Compute "on-site" densities (n1, ntild1, nhat1) for electron and positron ===== 3440 lmselect(:)=.true. 3441 opt_dens=1;if (include_nhat_in_gamma) opt_dens=0 3442 call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1,nspden_ep,1,& 3443 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij(iatom),& 3444 & pawtab(itypat),rho1,trho1) 3445 lmselect_ep(:)=.true. 3446 call pawdensities(rdum,cplex,iatom,lmselect_ep,lmselect_dum,lm_size,nhat1_ep,nspden_ep,1,& 3447 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_ep_(iatom),& 3448 & pawtab(itypat),rho1_ep,trho1_ep) 3449 !Compute contribution to annihilation rate 3450 3451 ABI_MALLOC(rhocore,(mesh_size)) 3452 3453 !First formalism: use densities on r,theta,phi 3454 if (dtset%pawxcdev==0) then 3455 3456 ABI_MALLOC(gamma,(mesh_size,2)) 3457 ABI_MALLOC(rhoarr1,(mesh_size)) 3458 ABI_MALLOC(rhoarr1_ep,(mesh_size)) 3459 3460 ! Loop on the angular part 3461 do ipt=1,pawang%angl_size 3462 ! Build densities 3463 rhoarr1=zero;rhoarr1_ep=zero;rhocore=zero 3464 do ilm=1,lm_size 3465 if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+rho1(:,ilm,1)*pawang%ylmr(ilm,ipt) 3466 end do 3467 do ilm=1,lm_size 3468 if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+rho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt) 3469 end do 3470 rhocore(:)=pawtab(itypat)%coredens(:) 3471 ! Make the densities positive 3472 if (electronpositron%particle==EP_ELECTRON) then 3473 call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos) 3474 call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos) 3475 else if (electronpositron%particle==EP_POSITRON) then 3476 call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1 ,dtset%xc_denpos) 3477 call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos) 3478 end if 3479 ! Compute Gamma 3480 ABI_MALLOC(grhoe2,(ngr)) 3481 ABI_MALLOC(grhocore2,(ngr)) 3482 if (electronpositron%particle==EP_ELECTRON) then 3483 call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,mesh_size,& 3484 & rhocore,rhoarr1_ep,rhoarr1,1) 3485 else if (electronpositron%particle==EP_POSITRON) then 3486 call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,mesh_size,& 3487 & rhocore,rhoarr1,rhoarr1_ep,1) 3488 end if 3489 ABI_FREE(grhoe2) 3490 ABI_FREE(grhocore2) 3491 ! Compute contribution to annihilation rates 3492 3493 ABI_MALLOC(ff,(mesh_size)) 3494 ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocorej(1:mesh_size) & 3495 & *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 3496 call simp_gen(intg,ff,pawrad(itypat)) 3497 intg=intg*pawang%angwgth(ipt)*four_pi 3498 rate =rate +intg 3499 ABI_FREE(ff) 3500 end do ! ipt 3501 ABI_FREE(gamma) 3502 ABI_FREE(rhoarr1) 3503 ABI_FREE(rhoarr1_ep) 3504 3505 !Second formalism: use (l,m) moments for densities 3506 else if (dtset%pawxcdev/=0) then 3507 3508 ! Build densities 3509 ABI_MALLOC(gammam,(mesh_size,lm_size)) 3510 ABI_MALLOC(rhotot,(mesh_size,lm_size)) 3511 ABI_MALLOC(rhotot_ep,(mesh_size,lm_size)) 3512 ABI_MALLOC(rhosph,(mesh_size)) 3513 ABI_MALLOC(rhosph_ep,(mesh_size)) 3514 3515 rhotot (:,:)=rho1 (:,:,1) 3516 rhotot_ep(:,:)=rho1_ep(:,:,1) 3517 rhocore(:)=pawtab(itypat)%coredens(:) 3518 rhosph (:)=rhotot (:,1)/sqfpi 3519 rhosph_ep(:)=rhotot_ep(:,1)/sqfpi 3520 ! Make spherical densities positive 3521 if (electronpositron%particle==EP_ELECTRON) then 3522 call mkdenpos(iwarnp,mesh_size,1,1,rhosph ,dtset%xc_denpos) 3523 call mkdenpos(iwarn ,mesh_size,1,1,rhosph_ep,dtset%xc_denpos) 3524 else if (electronpositron%particle==EP_POSITRON) then 3525 call mkdenpos(iwarn ,mesh_size,1,1,rhosph ,dtset%xc_denpos) 3526 call mkdenpos(iwarnp,mesh_size,1,1,rhosph_ep,dtset%xc_denpos) 3527 end if 3528 3529 ! Need gradients of electronic densities for GGA 3530 ABI_MALLOC(grhoe2,(ngr)) 3531 ABI_MALLOC(grhocore2,(ngr)) 3532 if (ngr>0) then 3533 if (electronpositron%particle==EP_ELECTRON) then 3534 call nderiv_gen(grhoe2,rhosph_ep,pawrad(itypat)) 3535 else if (electronpositron%particle==EP_POSITRON) then 3536 call nderiv_gen(grhoe2,rhosph,pawrad(itypat)) 3537 end if 3538 grhoe2(:)=grhoe2(:)**2 3539 call nderiv_gen(grhocore2,rhocore,pawrad(itypat)) 3540 grhocore2(:)=grhocore2(:)**2 3541 end if 3542 ! Compute Gamma for (rho-,rho+), 3543 ! (rho- +drho-,rho+), (rho- -drho-,rho+), 3544 ! (rho-,rho+ +drho+), (rho-,rho+ -drho+), 3545 ! (rho- +drho-,rho+ +drho+), (rho- -drho-,rho+ -drho+) 3546 ! Do a seven steps loop 3547 ABI_MALLOC(gam_,(mesh_size,2,7)) 3548 ABI_MALLOC(rho_,(mesh_size)) 3549 ABI_MALLOC(rho_ep_,(mesh_size)) 3550 ABI_MALLOC(rhocor_,(mesh_size)) 3551 ABI_MALLOC(grho2_,(ngr)) 3552 ABI_MALLOC(grhocor2_,(ngr)) 3553 3554 do ii=1,7 3555 ! Apply delta to get perturbed densities 3556 rho_(:)=rhosph(:);rho_ep_(:)=rhosph_ep(:);rhocor_(:)=rhocore(:) 3557 if (ngr>0) grho2_(:)=grhoe2(:) 3558 if (ngr>0) grhocor2_(:)=grhocore2(:) 3559 if (ii==2.or.ii==4.or.ii==6) fact=(one+delta) 3560 if (ii==3.or.ii==5.or.ii==7) fact=(one-delta) 3561 fact2=fact**2 3562 if (ii==2.or.ii==3.or.ii==6.or.ii==7) then 3563 rho_(:)=fact*rho_(:) 3564 if (electronpositron%particle==EP_POSITRON) then 3565 if (ngr>0) grho2_(:)=fact2*grho2_(:) 3566 rhocor_(:)=fact*rhocor_(:) 3567 if (ngr>0) grhocor2_(:)=fact2*grhocor2_(:) 3568 end if 3569 end if 3570 3571 if (ii==4.or.ii==5.or.ii==6.or.ii==7) then 3572 rho_ep_(:)=fact*rho_ep_(:) 3573 if (electronpositron%particle==EP_ELECTRON) then 3574 if (ngr>0) grho2_(:)=fact2*grho2_(:) 3575 rhocor_(:)=fact*rhocor_(:) 3576 if (ngr>0) grhocor2_(:)=fact2*grhocor2_(:) 3577 end if 3578 end if 3579 ! Compute gamma for these perturbed densities 3580 if (electronpositron%particle==EP_ELECTRON) then 3581 call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma,ngr,mesh_size,rhocor_,rho_ep_,rho_,1) 3582 else if (electronpositron%particle==EP_POSITRON) then 3583 call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma,ngr,mesh_size,rhocor_,rho_,rho_ep_,1) 3584 end if 3585 3586 end do ! end loop ii=1,7 3587 3588 ABI_FREE(rhocor_) 3589 ABI_FREE(grho2_) 3590 ABI_FREE(grhocor2_) 3591 ABI_FREE(grhoe2) 3592 ABI_FREE(grhocore2) 3593 rho_ (:)=rhosph (:);if (electronpositron%particle==EP_POSITRON) rho_ (:)=rho_ (:)+rhocore(:) 3594 rho_ep_(:)=rhosph_ep(:);if (electronpositron%particle==EP_ELECTRON) rho_ep_(:)=rho_ep_(:)+rhocore(:) 3595 ! Compute numerical first and second derivatives of Gamma 3596 ! d1gam(1) = dgam/drho+ (particle=ELECTRON), dgam/drho- (particle=POSITRON) 3597 ! d1gam(2) = dgam/drho- (particle=ELECTRON), dgam/drho+ (particle=POSITRON) 3598 ABI_MALLOC(d1gam,(mesh_size,2)) 3599 d1gam(:,:)=zero 3600 do ir=1,mesh_size 3601 if (rho_ (ir)>tol14) d1gam(ir,1)=(gam_(ir,1,2)-gam_(ir,1,3))*half/(delta*rho_ (ir)) 3602 if (rho_ep_ (ir)>tol14) d1gam(ir,2)=(gam_(ir,1,4)-gam_(ir,1,5))*half/(delta*rho_ep_ (ir)) 3603 end do 3604 3605 ! d2gam(1) = d2gam/drho+_drho+ (particle=ELECTRON), dgam/drho-_drho- (particle=POSITRON) 3606 ! d2gam(2) = d2gam/drho-_drho+ (particle=ELECTRON), dgam/drho+_drho- (particle=POSITRON) 3607 ! d2gam(3) = d2gam/drho-_drho- (particle=ELECTRON), dgam/drho+_drho+ (particle=POSITRON) 3608 ABI_MALLOC(d2gam,(mesh_size,3)) 3609 d2gam(:,:)=zero 3610 do ir=1,mesh_size 3611 if (rho_ (ir)>tol14) d2gam(ir,1)=(gam_(ir,1,2)+gam_(ir,1,3)-two*gam_(ir,1,1))/(delta*rho_ (ir))**2 3612 if (rho_ep_(ir)>tol14) then 3613 d2gam(ir,3)=(gam_(ir,1,4)+gam_(ir,1,5)-two*gam_(ir,1,1))/(delta*rho_ep_(ir))**2 3614 if (rho_(ir)>tol14) then 3615 d2gam(ir,2)=(gam_(ir,1,6)+gam_(ir,1,7)+two*gam_(ir,1,1) & 3616 & -gam_(ir,1,2)-gam_(ir,1,3)-gam_(ir,1,4)-gam_(ir,1,5)) & 3617 & *half/(delta*rho_(ir))/(delta*rho_ep_(ir)) 3618 end if 3619 end if 3620 end do 3621 3622 ABI_FREE(rho_) 3623 ABI_FREE(rho_ep_) 3624 ! Compute useful sums of densities 3625 ABI_MALLOC(v1sum,(mesh_size,3)) 3626 if ( dtset%pawxcdev>=2) then 3627 ABI_MALLOC(v2sum,(mesh_size,lm_size,3)) 3628 else 3629 ABI_MALLOC(v2sum,(0,0,0)) 3630 end if 3631 rhotot(:,1)=sqfpi*rhosph(:);rhotot_ep(:,1)=sqfpi*rhosph_ep(:) 3632 call pawxcsum(1,1,1,lmselect,lmselect_ep,lm_size,mesh_size,3,dtset%pawxcdev,& 3633 & pawang,rhotot,rhotot_ep,v1sum,v2sum) 3634 ! Compute final developpment of gamma moments 3635 gammam(:,:)=zero 3636 gammam(:,1)=gam_(:,1,1)*sqfpi 3637 gammam(:,1)=gammam(:,1)+(d2gam(:,2)*v1sum(:,2) & 3638 & +half*(d2gam(:,1)*v1sum(:,1)+d2gam(:,3)*v1sum(:,3)))/sqfpi 3639 do ilm=2,lm_size 3640 if (lmselect(ilm)) then 3641 gammam(:,ilm)=gammam(:,ilm)+d1gam(:,1)*rhotot(:,ilm) 3642 end if 3643 if (lmselect_ep(ilm)) then 3644 gammam(:,ilm)=gammam(:,ilm)+d1gam(:,2)*rhotot_ep(:,ilm) 3645 end if 3646 end do 3647 if (dtset%pawxcdev>1) then 3648 do ilm=2,lm_size 3649 gammam(:,ilm)=gammam(:,ilm)+d2gam(:,2)*v2sum(:,ilm,2) & 3650 & +half*(d2gam(:,1)*v2sum(:,ilm,1)+d2gam(:,3)*v2sum(:,ilm,3)) 3651 end do 3652 end if 3653 3654 ABI_FREE(gam_) 3655 ABI_FREE(d1gam) 3656 ABI_FREE(d2gam) 3657 ABI_FREE(v1sum) 3658 ABI_FREE(v2sum) 3659 ! Compute contribution to annihilation rate 3660 ABI_MALLOC(gg,(mesh_size,4)) 3661 gg=zero 3662 ABI_MALLOC(rhoarr2,(mesh_size)) 3663 do ilm=1,lm_size 3664 if (lmselect_ep(ilm)) gg(:,1)=gg(:,1)+rhotot_ep(:,ilm)*rhocorej(:)*gammam(:,ilm) 3665 end do 3666 ABI_FREE(rhoarr2) 3667 gg(1:mesh_size,1)=gg(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2 3668 call simp_gen(intg,gg(:,1),pawrad(itypat)) 3669 rate =rate +intg 3670 ABI_FREE(gg) 3671 ABI_FREE(gammam) 3672 ABI_FREE(rhotot) 3673 ABI_FREE(rhotot_ep) 3674 ABI_FREE(rhosph) 3675 ABI_FREE(rhosph_ep) 3676 3677 end if ! dtset%pawxcdev 3678 ABI_FREE(rhocore) 3679 3680 ABI_FREE(rho1) 3681 ABI_FREE(trho1) 3682 ABI_FREE(rho1_ep) 3683 ABI_FREE(trho1_ep) 3684 ABI_FREE(lmselect) 3685 ABI_FREE(lmselect_ep) 3686 ABI_FREE(lmselect_dum) 3687 ABI_FREE(nhat1) 3688 ABI_FREE(nhat1_ep) 3689 3690 !Reduction in case of distribution over atomic sites 3691 if (mpi_enreg%nproc_atom>1) then 3692 mpibuf=rate 3693 call xmpi_sum(mpibuf,mpi_enreg%comm_atom,ierr) 3694 rate=mpibuf 3695 end if 3696 3697 DBG_EXIT("COLL") 3698 3699 end subroutine posratecore
ABINIT/setup_positron [ Functions ]
NAME
setup_positron
FUNCTION
Do various initializations for the positron lifetime calculation
INPUTS
atindx(natom)=index table for atoms (see gstate.f) atindx1(natom)=index table for atoms, inverse of atindx dtefield <type(efield_type)> = variables related to Berry phase dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset ecore=core psp energy (part of total energy) (hartree) etotal=current value of total energy fock <type(fock_type)>= quantities to calculate Fock exact exchange forces_needed=if >0 forces are needed gred(3,natom)=gradients wrt nuclear positions in reduced coordinates gprimd(3,3)=dimensional primitive translations for reciprocal space gmet(3,3)=reciprocal space metric grchempottn(3,natom)=d(E_chemical_potential)/d(xred) (hartree) grcondft(3,natom)=d(E_constrainedDFT)/d(xred) (hartree) grewtn(3,natom)=d(Ewald)/d(xred) (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree) gsqcut=cutoff value on G**2 for sphere inside fft box hdr <type(hdr_type)>=the header of wf, den and pot files ifirst_gs= 0 if we are in a single ground-state calculation or in the first ground-state calculation of a structural minimization/dynamics indsym(4,nsym,natom)=index showing transformation of atom labels under symmetry operations (computed in symatm) istep=index of the number of steps in the routine scfcv istep_mix=index of the number of steps for the SCF mixing (can be <istep) kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 maxfor=maximum absolute value of fcart (forces) 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 mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor n3xccc=dimension of the xccc3d array (0 or nfftf). nattyp(ntypat)= # atoms of each type. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc nhat(nfftf,nspden*usepaw)= -PAW only- compensation density nkxc=second dimension of the array kxc, see rhotoxc.f for a description npwarr(nkpt)=number of planewaves in basis and on boundary for each k nvresid(nfftf,nspden)=array for the residual of the density/potential optres=0 if the potential residual has to be used for forces corrections =1 if the density residual has to be used for forces corrections paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawfgr(my_natom*usepaw) <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information (fine FFT grid) ph1dc(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases (coarse FFT grid) psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translations in real space (bohr) stress_needed=if >0 stresses are needed strscondft(6)=cDFT correction to stress strsxc(6)=xc correction to stress symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates ucvol=unit cell volume in bohr**3. usecprj= 1 if cprj array is stored in memory vhartr(nfftf)=array for holding Hartree potential vpsp(nfftf)=array for holding local psp vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space vxctau(nfft,nspden,4*usekden)=(only for meta-GGA) derivative of XC energy density wrt kinetic energy density (depsxcdtau) xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xcctau3d(n3xccc*usekden)=(only for meta-GGA): 3D core electron kinetic energy density for XC core correction xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation energies <type(energies_type)>=all part of total energy. cg(2,mcg)=wavefunctions cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors: cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector. eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) occ(mband*nkpt*nsppol)=occupation number for each band at each k point pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data rhog(2,nfft)=Fourier transform of total electron/positron density rhor(nfft,nspden)=total electron/positron density (el/bohr**3)
SOURCE
181 subroutine setup_positron(atindx,atindx1,cg,cprj,dtefield,dtfil,dtset,ecore,eigen,etotal,electronpositron,& 182 & energies,fock,forces_needed,gred,gmet,gprimd,grchempottn,& 183 & grcondft,grewtn,grvdw,gsqcut,hdr,extfpmd,ifirst_gs,indsym,istep,istep_mix,kg,& 184 & kxc,maxfor,mcg,mcprj,mgfft,mpi_enreg,my_natom,n3xccc,nattyp,nfft,ngfft,ngrvdw,nhat,nkxc,npwarr,nvresid,occ,optres,& 185 & paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1dc,psps,rhog,rhor,& 186 & rmet,rprimd,stress_needed,strscondft,strsxc,symrec,ucvol,usecprj,vhartr,vpsp,vxc,vxctau,& 187 & xccc3d,xcctau3d,xred,ylm,ylmgr) 188 189 !Arguments ------------------------------------ 190 !scalars 191 integer,intent(in) :: forces_needed,ifirst_gs,istep,mcg,mcprj,mgfft,my_natom,n3xccc,nfft 192 integer,intent(in) :: ngrvdw,nkxc,optres,stress_needed,usecprj 193 integer,intent(inout) :: istep_mix 194 real(dp),intent(in) :: ecore,etotal,gsqcut,maxfor,ucvol 195 type(efield_type),intent(in) :: dtefield 196 type(datafiles_type),intent(in) :: dtfil 197 type(dataset_type),intent(in) :: dtset 198 type(electronpositron_type),pointer :: electronpositron 199 type(energies_type),intent(inout) :: energies 200 type(hdr_type),intent(inout) :: hdr 201 type(extfpmd_type),pointer,intent(inout) :: extfpmd 202 type(MPI_type),intent(inout) :: mpi_enreg 203 type(pawang_type),intent(in) :: pawang 204 type(pawfgr_type),intent(in) :: pawfgr 205 type(pseudopotential_type), intent(in) :: psps 206 type(fock_type),pointer, intent(inout) :: fock 207 !arrays 208 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom) 209 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%natom),ngfft(18) 210 integer,intent(in) :: npwarr(dtset%nkpt),symrec(3,3,dtset%nsym) 211 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),grchempottn(3,dtset%natom),grcondft(3,dtset%natom) 212 real(dp),intent(in) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc) 213 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),ph1dc(2,(3*(2*dtset%mgfft+1)*dtset%natom)*dtset%usepaw) 214 real(dp),intent(in) :: rmet(3,3),strscondft(6),strsxc(6),vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden) 215 real(dp),intent(in) :: vxctau(nfft,dtset%nspden,4*dtset%usekden) 216 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 217 real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 218 real(dp),intent(inout) :: cg(2,mcg) 219 real(dp),intent(inout) :: nhat(nfft,dtset%nspden*dtset%usepaw) 220 real(dp),intent(inout) :: nvresid(nfft,dtset%nspden) 221 real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol),gred(3,dtset%natom) 222 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 223 real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,dtset%nspden),rprimd(3,3) 224 real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*dtset%usekden),xred(3,dtset%natom) 225 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj) 226 type(paw_ij_type),intent(in) :: paw_ij(my_natom*dtset%usepaw) 227 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*dtset%usepaw) 228 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw) 229 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 230 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*dtset%usepaw) 231 232 !Local variables------------------------------- 233 !scalars 234 integer,parameter :: cplex1=1 235 integer :: history_level,iatom,iband,icalctype,icalctype0,icg,ifft,ikpt 236 integer :: iocc,ireadwf,ispden,isppol,n3xccc0,nocc,optfor,optstr,rdwrpaw,comm_cell 237 logical,parameter :: always_restart=.false. ! Set to true to restart by a pure electronic step at each new atomic structure 238 logical :: need_scocc,new_calctype 239 real(dp) :: boxcut_dum,diffor_dum,ecut_eff,eigtmp,etotal_read,gsqcut_eff,maxfor_dum 240 real(dp) :: maxocc,nelect,occlast,occtmp,rhotmp 241 character(len=69) :: TypeCalcStrg 242 character(len=500) :: message 243 character(len=fnlen) :: fname 244 type(energies_type) :: energies_tmp 245 type(wvl_data) :: wvl 246 type(hdr_type) :: hdr_den 247 !arrays 248 integer,allocatable :: nlmn(:) 249 real(dp) :: cgtmp(2) 250 real(dp),parameter :: qphon(3)=(/zero,zero,zero/) 251 real(dp),allocatable :: favg_dum(:),fcart_dum(:,:),forold_dum(:,:),gred_tmp(:,:) 252 real(dp),allocatable :: gresid_dum(:,:),grhf_dum(:,:),grxc_dum(:,:) 253 real(dp),allocatable :: rhog_ep(:,:),scocc(:),str_tmp(:),synlgr_dum(:,:) 254 real(dp) :: nhatgr(0,0,0) 255 type(pawcprj_type),allocatable :: cprj_tmp(:,:) 256 type(pawrhoij_type),allocatable :: pawrhoij_tmp(:) 257 258 ! ************************************************************************* 259 260 !Compatibility tests 261 if (dtset%positron==0) then 262 ABI_BUG('Not valid for dtset%positron=0!') 263 end if 264 265 if (istep>1.and.nfft/=electronpositron%nfft) then 266 ABI_BUG('Invalid value for nfft!') 267 end if 268 269 if (dtset%usewvl==1) then 270 ABI_BUG('Not valid for wavelets!') 271 end if 272 273 if (dtset%positron==1) then 274 do isppol=1,dtset%nsppol 275 do ikpt=1,dtset%nkpt 276 if (dtset%nband(ikpt+dtset%nkpt*(isppol-1))/=dtset%nband(1)) then 277 message = "dtset%positron needs nband to be the same at each k-point !" 278 ABI_ERROR(message) 279 end if 280 end do 281 end do 282 end if 283 284 comm_cell = mpi_enreg%comm_cell 285 286 !----------------------------------------------------------------------- 287 !Compute new value for calctype (type of electron-positron calculation) 288 !----------------------------------------------------------------------- 289 icalctype0=electronpositron%calctype 290 291 new_calctype=.false. 292 if (dtset%positron==1.or.dtset%positron==2) then 293 if (istep==1) new_calctype=.true. 294 else if (dtset%positron<0) then 295 if (ifirst_gs/=0.and.istep==1.and.(.not.always_restart)) new_calctype=.true. 296 if (electronpositron%scf_converged) new_calctype=.true. 297 end if 298 299 !Comment: 300 !history_level=-1: not used 301 !history_level= 0: rhor from scratch, rhor_ep from scratch or read 302 !history_level= 1: rhor in memory, rhor_ep from scratch or read 303 !history_level= 2: rhor_ep <-rhor, rhor from scratch 304 !history_level= 3: rhor_ep <-> rhor 305 !history_level= 4: rhor in memory, rhor_ep in memory 306 history_level=-1 307 if (dtset%positron==1.or.dtset%positron==2) then 308 if (ifirst_gs==0.and.istep==1) history_level=0 309 if (ifirst_gs/=0.and.istep==1) history_level=4 310 else if (dtset%positron<0) then 311 if (.not.electronpositron%scf_converged) then 312 if (ifirst_gs/=0.and.istep==1.and.(.not.always_restart)) history_level=4 313 else if (electronpositron%scf_converged) then 314 if (icalctype0==0) history_level=2 315 if (icalctype0> 0) history_level=3 316 end if 317 end if 318 319 electronpositron%calctype=icalctype0 320 if (dtset%positron==1.or.dtset%positron==2) then 321 electronpositron%calctype=dtset%positron 322 else if (dtset%positron<0) then 323 if (electronpositron%scf_converged) then 324 if (icalctype0==0) electronpositron%calctype=1 325 if (icalctype0>0 ) electronpositron%calctype=3-electronpositron%calctype 326 else if (ifirst_gs/=0.and.istep==1) then 327 if (always_restart) then 328 electronpositron%calctype=0 329 else 330 electronpositron%calctype=2 331 ! if (electronpositron%particle==EP_POSITRON) electronpositron%calctype=1 332 end if 333 end if 334 end if 335 336 electronpositron%scf_converged=.false. 337 if (new_calctype) electronpositron%istep=electronpositron%istep+1 338 if (istep==1) electronpositron%istep=1 339 ireadwf=dtfil%ireadwf;if (electronpositron%istep>1) ireadwf=0 340 341 !============================================ 342 !The following lines occur only when the type 343 !of electron-positron calculation changes 344 !============================================ 345 if (new_calctype) then 346 347 ! Reset some indexes 348 if (electronpositron%calctype==0) then 349 electronpositron%particle=EP_NOTHING 350 else if (electronpositron%calctype==1) then 351 electronpositron%particle=EP_ELECTRON 352 else if (electronpositron%calctype==2) then 353 electronpositron%particle=EP_POSITRON 354 end if 355 electronpositron%has_pos_ham=mod(electronpositron%calctype,2) 356 electronpositron%istep_scf=1 357 istep_mix=1 358 359 ! ----------------------------------------------------------------------------------------- 360 ! Update forces and stresses 361 ! If electronpositron%calctype==1: gred_ep/stress_ep are the electronic gred/stress 362 ! If electronpositron%calctype==2: gred_ep/stress_ep are the positronic gred/stress 363 ! ----------------------------------------------------------------------------------------- 364 if (history_level==2.or.history_level==3) then 365 optstr=0;optfor=0 366 if (allocated(electronpositron%stress_ep)) optstr=stress_needed 367 if (allocated(electronpositron%gred_ep).and.forces_needed==2) optfor=1 368 if (optfor>0.or.optstr>0) then 369 ABI_MALLOC(favg_dum,(3)) 370 ABI_MALLOC(fcart_dum,(3,dtset%natom)) 371 ABI_MALLOC(forold_dum,(3,dtset%natom)) 372 ABI_MALLOC(gresid_dum,(3,dtset%natom)) 373 ABI_MALLOC(grhf_dum,(3,dtset%natom)) 374 ABI_MALLOC(grxc_dum,(3,dtset%natom)) 375 ABI_MALLOC(synlgr_dum,(3,dtset%natom)) 376 ABI_MALLOC(gred_tmp,(3,dtset%natom)) 377 ABI_MALLOC(str_tmp,(6)) 378 forold_dum=zero;n3xccc0=n3xccc 379 icalctype=electronpositron%calctype;electronpositron%calctype=-icalctype0 380 if (electronpositron%calctype==0) electronpositron%calctype=-100 381 if (electronpositron%calctype==-1) n3xccc0=0 ! Note: if calctype=-1, previous calculation was positron 382 call forstr(atindx1,cg,cprj,diffor_dum,dtefield,dtset,eigen,electronpositron,energies,& 383 & favg_dum,fcart_dum,fock,forold_dum,gred_tmp,grchempottn,grcondft,gresid_dum,grewtn,grhf_dum,grvdw,grxc_dum,gsqcut,& 384 & extfpmd,indsym,kg,kxc,maxfor_dum,mcg,mcprj,mgfft,mpi_enreg,my_natom,n3xccc0,nattyp,nfft,ngfft,& 385 & ngrvdw,nhat,nkxc,npwarr,dtset%ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,& 386 & pawfgrtab,pawrad,pawrhoij,pawtab,ph1dc,ph1d,psps,rhog,rhor,rprimd,optstr,strscondft,strsxc,str_tmp,symrec,& 387 & synlgr_dum,ucvol,usecprj,vhartr,vpsp,vxc,vxctau,wvl,xccc3d,xcctau3d,xred,ylm,ylmgr,0.0_dp) 388 electronpositron%calctype=icalctype 389 if (optfor>0) electronpositron%gred_ep(:,:)=gred_tmp(:,:) 390 if (optstr>0) electronpositron%stress_ep(:)=str_tmp(:) 391 ABI_FREE(favg_dum) 392 ABI_FREE(fcart_dum) 393 ABI_FREE(forold_dum) 394 ABI_FREE(gresid_dum) 395 ABI_FREE(grhf_dum) 396 ABI_FREE(grxc_dum) 397 ABI_FREE(synlgr_dum) 398 ABI_FREE(gred_tmp) 399 ABI_FREE(str_tmp) 400 end if 401 if (optfor==0.and.forces_needed>0.and.allocated(electronpositron%gred_ep)) then 402 electronpositron%gred_ep(:,:)=gred(:,:)-electronpositron%gred_ep(:,:) 403 end if 404 end if 405 406 ! ---------------------------------------------------------------------------------------------------- 407 ! Initialize/Update densities 408 ! If electronpositron%calctype==1: rhor is the positronic density, rhor_ep is the electronic density 409 ! If electronpositron%calctype==2: rhor is the electronic density, rhor_ep is the positronic density 410 ! --------------------------------------------------------------------------------------------------- 411 ABI_MALLOC(rhog_ep,(2,nfft)) 412 413 ! ===== PREVIOUS DENSITY RHOR_EP: 414 if (history_level==0.or.history_level==1) then 415 ! ----- Read from disk 416 if (dtset%positron>0) then 417 rdwrpaw=dtset%usepaw 418 fname=trim(dtfil%fildensin);if (dtset%positron==2) fname=trim(dtfil%fildensin)//'_POSITRON' 419 call read_rhor(trim(fname), cplex1, dtset%nspden, nfft, ngfft, rdwrpaw, mpi_enreg, electronpositron%rhor_ep, & 420 hdr_den, electronpositron%pawrhoij_ep, comm_cell, check_hdr=hdr) 421 etotal_read = hdr_den%etot; call hdr_den%free() 422 call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,1,ngfft,0) 423 if (dtset%usepaw==1.and.allocated(electronpositron%nhat_ep)) then 424 call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,& 425 & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,electronpositron%nhat_ep,& 426 & electronpositron%pawrhoij_ep,electronpositron%pawrhoij_ep,pawtab,& 427 & qphon,rprimd,ucvol,dtset%usewvl,xred,distribfft=mpi_enreg%distribfft,& 428 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 429 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0) 430 end if 431 end if 432 ! ----- Electronic from scratch 433 if (dtset%positron<0.and.electronpositron%calctype==1) then 434 ecut_eff=dtset%pawecutdg*(dtset%dilatmx)**2 435 call getcut(boxcut_dum,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,std_out,qphon,ngfft) 436 call initro(atindx,dtset%densty,gmet,gsqcut_eff,dtset%usepaw,mgfft,mpi_enreg,& 437 & psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,dtset%ntypat,& 438 & psps,pawtab,ph1d,psps%qgrid_vl,rhog_ep,electronpositron%rhor_ep,& 439 & dtset%spinat,ucvol,dtset%usepaw,dtset%ziontypat,dtset%znucl) 440 if (dtset%usepaw==1) then 441 if (size(electronpositron%pawrhoij_ep)>0) then 442 ABI_MALLOC(pawrhoij_tmp,(my_natom)) 443 call initrhoij(electronpositron%pawrhoij_ep(1)%cplex_rhoij,dtset%lexexch,& 444 & dtset%lpawu,my_natom,dtset%natom,dtset%nspden,& 445 & electronpositron%pawrhoij_ep(1)%nspinor,dtset%nsppol,& 446 & dtset%ntypat,pawrhoij_tmp,dtset%pawspnorb,pawtab,cplex1,dtset%spinat,dtset%typat,& 447 & ngrhoij=electronpositron%pawrhoij_ep(1)%ngrhoij,& 448 & nlmnmix=electronpositron%pawrhoij_ep(1)%lmnmix_sz,& 449 & use_rhoij_=electronpositron%pawrhoij_ep(1)%use_rhoij_,& 450 & use_rhoijres=electronpositron%pawrhoij_ep(1)%use_rhoijres,& 451 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 452 if (electronpositron%pawrhoij_ep(1)%lmnmix_sz>0) then 453 do iatom=1,my_natom 454 pawrhoij_tmp(iatom)%kpawmix(:)=electronpositron%pawrhoij_ep(iatom)%kpawmix(:) 455 end do 456 end if 457 call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep) 458 call pawrhoij_free(pawrhoij_tmp) 459 ABI_FREE(pawrhoij_tmp) 460 end if 461 if (allocated(electronpositron%nhat_ep)) then 462 call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,& 463 & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,electronpositron%nhat_ep,& 464 & electronpositron%pawrhoij_ep,electronpositron%pawrhoij_ep,pawtab,qphon,rprimd,& 465 & ucvol,dtset%usewvl,xred,distribfft=mpi_enreg%distribfft,& 466 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 467 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0) 468 end if 469 end if 470 end if 471 ! ----- Positronic from scratch 472 if (dtset%positron<0.and.electronpositron%calctype==2) then 473 electronpositron%rhor_ep(:,1)=one/ucvol 474 if (dtset%nspden>=2) electronpositron%rhor_ep(:,2)=half/ucvol 475 if (dtset%nspden==4) electronpositron%rhor_ep(:,3:4)=zero 476 rhog_ep=zero;rhog_ep(1,1)=one/ucvol 477 if (dtset%usepaw==1) then 478 do iatom=1,dtset%natom 479 electronpositron%pawrhoij_ep(iatom)%rhoijp(:,:)=zero 480 electronpositron%pawrhoij_ep(iatom)%nrhoijsel=0 481 end do 482 if (allocated(electronpositron%nhat_ep)) electronpositron%nhat_ep=zero 483 end if 484 end if 485 end if 486 ! ----- Deduced from rhor in memory 487 if (history_level==2) then 488 electronpositron%rhor_ep(:,:)=rhor(:,:) 489 rhog_ep(:,:)=rhog(:,:) 490 if (dtset%usepaw==1) then 491 call pawrhoij_copy(pawrhoij,electronpositron%pawrhoij_ep) 492 if (allocated(electronpositron%nhat_ep)) electronpositron%nhat_ep(:,:)=nhat(:,:) 493 end if 494 end if 495 496 ! ===== CURRENT DENSITY RHOR: 497 if (history_level==0.or.history_level==2) then 498 if (ireadwf==0) then 499 ! ----- Positronic from scratch 500 if (electronpositron%calctype==1) then 501 rhor(:,1)=one/ucvol 502 if (dtset%nspden>=2) rhor(:,2)=half/ucvol 503 if (dtset%nspden==4) rhor(:,3:4)=zero 504 rhog=zero;rhog(1,1)=one/ucvol 505 if (dtset%usepaw==1) then 506 do iatom=1,my_natom 507 pawrhoij(iatom)%rhoijp(:,:)=zero 508 pawrhoij(iatom)%nrhoijsel=0 509 end do 510 nhat(:,:)=zero 511 end if 512 end if 513 ! ----- Electronic from scratch 514 if (electronpositron%calctype==2) then 515 ecut_eff=dtset%pawecutdg*(dtset%dilatmx)**2 516 call getcut(boxcut_dum,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,std_out,qphon,ngfft) 517 call initro(atindx,dtset%densty,gmet,gsqcut_eff,dtset%usepaw,mgfft,mpi_enreg,& 518 & psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,dtset%ntypat,& 519 & psps,pawtab,ph1d,psps%qgrid_vl,rhog,rhor,dtset%spinat,ucvol,& 520 & dtset%usepaw,dtset%ziontypat,dtset%znucl) 521 522 if (dtset%usepaw==1) then 523 if (size(pawrhoij)>0) then 524 ABI_MALLOC(pawrhoij_tmp,(my_natom)) 525 call initrhoij(pawrhoij(1)%cplex_rhoij,dtset%lexexch,dtset%lpawu,& 526 & my_natom,dtset%natom,dtset%nspden,pawrhoij(1)%nspinor,dtset%nsppol,& 527 & dtset%ntypat,pawrhoij_tmp,dtset%pawspnorb,pawtab,pawrhoij(1)%qphase,dtset%spinat,& 528 & dtset%typat,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,& 529 & use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres,& 530 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 531 do iatom=1,my_natom 532 pawrhoij_tmp(iatom)%kpawmix(:)=pawrhoij(iatom)%kpawmix(:) 533 end do 534 call pawrhoij_copy(pawrhoij_tmp,pawrhoij) 535 call pawrhoij_free(pawrhoij_tmp) 536 ABI_FREE(pawrhoij_tmp) 537 end if 538 call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,& 539 & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,nhat,& 540 & pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, & 541 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 542 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,& 543 & me_g0=mpi_enreg%me_g0,distribfft=mpi_enreg%distribfft) 544 end if 545 546 end if 547 end if 548 end if 549 550 ! ===== EXCHANGE POSITRONIC AND ELECTRONIC DENSITY (CURRENT AND PREVIOUS) 551 if (history_level==3) then 552 do ispden=1,dtset%nspden 553 do ifft=1,nfft 554 rhotmp=rhor(ifft,ispden) 555 rhor(ifft,ispden)=electronpositron%rhor_ep(ifft,ispden) 556 electronpositron%rhor_ep(ifft,ispden)=rhotmp 557 end do 558 end do 559 rhog_ep(:,:)=rhog 560 call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,1,ngfft,0) 561 ! If PAW, exchange "positronic" and "electronic" rhoij 562 if (dtset%usepaw==1) then 563 if (size(pawrhoij)>0.and.size(electronpositron%pawrhoij_ep)>0) then 564 ABI_MALLOC(pawrhoij_tmp,(my_natom)) 565 call pawrhoij_alloc(pawrhoij_tmp,pawrhoij(1)%cplex_rhoij,pawrhoij(1)%nspden,& 566 & pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,dtset%typat,& 567 & pawtab=pawtab,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,& 568 & use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres, & 569 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 570 call pawrhoij_copy(pawrhoij,pawrhoij_tmp) 571 call pawrhoij_copy(electronpositron%pawrhoij_ep,pawrhoij) 572 call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep) 573 call pawrhoij_free(pawrhoij_tmp) 574 ABI_FREE(pawrhoij_tmp) 575 end if 576 if (allocated(electronpositron%nhat_ep)) then 577 do ispden=1,dtset%nspden 578 do ifft=1,nfft 579 rhotmp=nhat(ifft,ispden) 580 nhat(ifft,ispden)=electronpositron%nhat_ep(ifft,ispden) 581 electronpositron%nhat_ep(ifft,ispden)=rhotmp 582 end do 583 end do 584 end if 585 end if 586 end if 587 588 ! ===== COMPUTE HARTREE POTENTIAL ASSOCIATED TO RHOR_EP 589 if (history_level==4) then 590 call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,1,ngfft,0) 591 end if 592 if (history_level/=-1) then 593 call hartre(1,gsqcut,dtset%icutcoul,dtset%usepaw,mpi_enreg,nfft,ngfft,& 594 & dtset%nkpt,dtset%rcut,rhog_ep,rprimd,dtset%vcutgeo,electronpositron%vha_ep) 595 electronpositron%vha_ep=-electronpositron%vha_ep 596 else 597 electronpositron%vha_ep=zero 598 end if 599 ABI_FREE(rhog_ep) 600 601 ! ---------------------------------------------------------------------- 602 ! Initialize/Update energies 603 ! ---------------------------------------------------------------------- 604 electronpositron%etotal_prev=etotal 605 electronpositron%maxfor_prev=maxfor 606 607 ! Inits/exchange news energies 608 ! Retrieve energy of non-evolving particle(s) 609 if (history_level== 0) then 610 call energies_init(energies) 611 call energies_init(electronpositron%energies_ep) 612 if (dtset%positron>0) energies%e0_electronpositron=etotal_read 613 if (dtset%positron<0) energies%e0_electronpositron=zero 614 else if (history_level== 1) then 615 call energies_init(electronpositron%energies_ep) 616 if (dtset%positron>0) energies%e0_electronpositron=etotal_read 617 else if (history_level== 2) then 618 call energies_copy(energies,electronpositron%energies_ep) 619 call energies_init(energies) 620 energies%e0_electronpositron=electronpositron%e0 621 else if (history_level== 3) then 622 call energies_copy(electronpositron%energies_ep,energies_tmp) 623 call energies_copy(energies,electronpositron%energies_ep) 624 call energies_copy(energies_tmp,energies) 625 energies%e0_electronpositron=electronpositron%e0 626 ! else if (history_level== 4) then 627 end if 628 629 ! Adjust core psps energy 630 if (electronpositron%calctype/=1) energies%e_corepsp=ecore/ucvol 631 632 ! ----------------------------------------------------------------------------------------- 633 ! Update wavefunctions 634 ! If electronpositron%calctype==1: cg are the positronic WFs, cg_ep are the electronic WFs 635 ! If electronpositron%calctype==2: cg are the electronic WFs, cg_ep are the positronic WFs 636 ! ----------------------------------------------------------------------------------------- 637 if (electronpositron%dimcg>0.or.electronpositron%dimcprj>0) then 638 639 if (history_level==0.or.history_level==1) then 640 electronpositron%cg_ep=zero 641 end if 642 643 if (history_level==2) then 644 if (electronpositron%dimcg>0) then 645 do icg=1,electronpositron%dimcg 646 electronpositron%cg_ep(1:2,icg)=cg(1:2,icg) 647 end do 648 end if 649 if (dtset%usepaw==1.and.electronpositron%dimcprj>0) then 650 call pawcprj_copy(cprj,electronpositron%cprj_ep) 651 end if 652 end if 653 654 if (history_level==3) then 655 if (electronpositron%dimcg>0) then 656 do icg=1,electronpositron%dimcg 657 cgtmp(1:2)=electronpositron%cg_ep(1:2,icg) 658 electronpositron%cg_ep(1:2,icg)=cg(1:2,icg) 659 cg(1:2,icg)=cgtmp(1:2) 660 end do 661 end if 662 if (dtset%usepaw==1.and.electronpositron%dimcprj>0) then 663 ABI_MALLOC(nlmn,(dtset%natom)) 664 ABI_MALLOC(cprj_tmp,(dtset%natom,electronpositron%dimcprj)) 665 do iatom=1,dtset%natom 666 nlmn(iatom)=cprj(iatom,1)%nlmn 667 end do 668 call pawcprj_alloc(cprj_tmp,cprj(1,1)%ncpgr,nlmn) 669 ABI_FREE(nlmn) 670 call pawcprj_copy(electronpositron%cprj_ep,cprj_tmp) 671 call pawcprj_copy(cprj,electronpositron%cprj_ep) 672 call pawcprj_copy(cprj_tmp,cprj) 673 call pawcprj_free(cprj_tmp) 674 ABI_FREE(cprj_tmp) 675 end if 676 end if 677 678 end if ! dimcg>0 or dimcprj>0 679 680 ! ----------------------------------------------------------------------------------------------------------- 681 ! Initialize/Update occupations 682 ! If electronpositron%calctype==1: occ are the positronic occupations, occ_ep are the electronic occupations 683 ! If electronpositron%calctype==2: occ are the electronic occupations, occ_ep are the positronic occupations 684 ! ----------------------------------------------------------------------------------------------------------- 685 ! When needed, precompute electronic occupations with semiconductor occupancies 686 need_scocc=.false. 687 if (electronpositron%dimocc>0.and.electronpositron%calctype==1.and. & 688 & (history_level==0.or.history_level==1)) need_scocc=.true. 689 if (electronpositron%calctype==2.and.ireadwf==0.and. & 690 & (history_level==0.or.history_level==2.or. & 691 & (history_level==3.and.electronpositron%dimocc==0))) need_scocc=.true. 692 if (need_scocc) then 693 nelect=-dtset%cellcharge(1) 694 do iatom=1,dtset%natom 695 nelect=nelect+dtset%ziontypat(dtset%typat(iatom)) 696 end do 697 maxocc=two/real(dtset%nsppol*dtset%nspinor,dp) 698 nocc=int((nelect-tol8)/maxocc) + 1 699 nocc=min(nocc,dtset%nband(1)*dtset%nsppol) 700 occlast=nelect-maxocc*(nocc-1) 701 ABI_MALLOC(scocc,(dtset%nband(1)*dtset%nsppol)) 702 scocc=zero 703 if (1<nocc) scocc(1:nocc-1)=maxocc 704 if (1<=nocc) scocc(nocc)=occlast 705 end if 706 707 ! ===== PREVIOUS OCCUPATIONS OCC_EP: 708 if (electronpositron%dimocc>0) then 709 if (history_level==0.or.history_level==1) then 710 ! ----- Electronic from scratch 711 if (electronpositron%calctype==1) then 712 ! Initialize electronic occupations with semiconductor occupancies 713 do ikpt=1,dtset%nkpt 714 do iband=1,dtset%nband(1) 715 do isppol=1,dtset%nsppol 716 electronpositron%occ_ep(iband+dtset%nband(1)*(ikpt-1+dtset%nkpt*(isppol-1)))=& 717 & scocc(isppol+dtset%nsppol*(iband-1)) 718 end do 719 end do 720 end do 721 end if 722 ! ----- Positronic from scratch 723 if (electronpositron%calctype==1) then 724 ! Initialize positronic occupations with only one positron (or less) 725 electronpositron%occ_ep(:)=zero 726 isppol=1;iocc=1 727 do ikpt=1,dtset%nkpt 728 electronpositron%occ_ep(iocc)=electronpositron%posocc 729 iocc=iocc+dtset%nband(ikpt+dtset%nkpt*(isppol-1)) 730 end do 731 end if 732 end if 733 ! ----- Deduced from occ in memory 734 if (history_level==2) then 735 electronpositron%occ_ep(:)=occ(:) 736 end if 737 end if ! dimocc>0 738 739 ! ===== CURRENT OCCUPATIONS OCC: 740 if (history_level==0.or.history_level==2.or.(history_level==3.and.electronpositron%dimocc==0)) then 741 if (ireadwf==0) then 742 ! ----- Positronic from scratch 743 if (electronpositron%calctype==1) then 744 ! Initialize positronic occupations with only one positron (or less) 745 occ(:)=zero 746 isppol=1;iocc=1 747 do ikpt=1,dtset%nkpt 748 occ(iocc)=electronpositron%posocc 749 iocc=iocc+dtset%nband(ikpt+dtset%nkpt*(isppol-1)) 750 end do 751 end if 752 ! ----- Electronic from scratch 753 if (electronpositron%calctype==2) then 754 ! Initialize electronic occupations with semiconductor occupancies 755 do ikpt=1,dtset%nkpt 756 do iband=1,dtset%nband(1) 757 do isppol=1,dtset%nsppol 758 occ(iband+dtset%nband(1)*(ikpt-1+dtset%nkpt*(isppol-1)))=& 759 & scocc(isppol+dtset%nsppol*(iband-1)) 760 end do 761 end do 762 end do 763 end if 764 end if 765 end if 766 767 ! ===== EXCHANGE POSITRONIC AND ELECTRONIC OCCUPATIONS (CURRENT AND PREVIOUS) 768 if (history_level==3.and.electronpositron%dimocc>0) then 769 do iocc=1,electronpositron%dimocc 770 occtmp=occ(iocc) 771 occ(iocc)=electronpositron%occ_ep(iocc) 772 electronpositron%occ_ep(iocc)=occtmp 773 end do 774 end if 775 776 if (need_scocc) then 777 ABI_FREE(scocc) 778 end if 779 780 ! ----------------------------------------------------------------------------------------------------------- 781 ! Initialize/Update eigen energies 782 ! If electronpositron%calctype==1: eigen are the positronic eigen E, eigen_ep are the electronic eigen E 783 ! If electronpositron%calctype==2: eigen are the electronic eigen E, eigen_ep are the positronic eigen E 784 ! ----------------------------------------------------------------------------------------------------------- 785 786 ! ===== PREVIOUS EIGEN ENERGIES EIGEN_EP: 787 if (electronpositron%dimeigen>0) then 788 if (history_level==0.or.history_level==1) then 789 ! ----- Electronic or positronic from scratch 790 electronpositron%eigen_ep(:)=zero 791 end if 792 ! ----- Deduced from eigen in memory 793 if (history_level==2) then 794 electronpositron%eigen_ep(:)=eigen(:) 795 end if 796 end if ! dimeigen>0 797 798 ! ===== CURRENT EIGEN ENERGIES EIGEN: 799 if (history_level==0.or.history_level==2.or.(history_level==3.and.electronpositron%dimeigen==0)) then 800 if (ireadwf==0) then 801 ! ----- Electronic or positronic from scratch 802 eigen(:)=zero 803 end if 804 end if 805 806 ! ===== EXCHANGE POSITRONIC AND ELECTRONIC EIGEN ENERGIES (CURRENT AND PREVIOUS) 807 if (history_level==3.and.electronpositron%dimeigen>0) then 808 do iocc=1,electronpositron%dimeigen 809 eigtmp=eigen(iocc) 810 eigen(iocc)=electronpositron%eigen_ep(iocc) 811 electronpositron%eigen_ep(iocc)=eigtmp 812 end do 813 end if 814 815 ! In some cases cprj are kept in memory, so we have to update them before the call of vtorho 816 if (dtset%cprj_in_memory/=0) then 817 iatom=0 818 call wrtout(std_out,' Computing cprj from wavefunctions (positron)') 819 call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,0,& 820 & 0,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,& 821 & dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,& 822 & dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,& 823 & xred,ylm,ylmgr) 824 call wrtout(std_out,' cprj is computed') 825 end if 826 827 ! ============================================= 828 end if ! the type of e-p calculation changes 829 !============================================= 830 831 !------------------------------------------------------------------ 832 !Write messages 833 !------------------------------------------------------------------ 834 if (istep_mix==1.and.dtset%positron/=0) then 835 ! Log message 836 if (electronpositron%calctype==0) then 837 message = 'Were are now performing an electronic ground-state calculation...' 838 else if (electronpositron%calctype==1) then 839 message = 'Were are now performing a positronic ground-state calculation...' 840 else if (electronpositron%calctype==2) then 841 message = 'Were are now performing an electronic ground-state calculation in presence of a positron...' 842 end if 843 ABI_COMMENT(message) 844 ! Output message 845 if (dtset%positron<0) then 846 if (electronpositron%calctype==0) then 847 TypeCalcStrg='ELECTRONIC GROUND-STATE CALCULATION' 848 else if (electronpositron%calctype==1) then 849 TypeCalcStrg='POSITRONIC GROUND-STATE CALCULATION IN PRESENCE OF ELECTRONS AND IONS' 850 else if (electronpositron%calctype==2) then 851 TypeCalcStrg='ELECTRONIC GROUND-STATE CALCULATION IN PRESENCE OF A POSITRON' 852 end if 853 if (istep>1) then 854 write(message,'(2a,i3,2a)') ch10,'TC-DFT STEP ',electronpositron%istep,' - ',trim(TypeCalcStrg) 855 else 856 write(message,'(a,i3,2a)') 'TC-DFT STEP ',electronpositron%istep,' - ',trim(TypeCalcStrg) 857 end if 858 call wrtout(ab_out,message,'COLL') 859 end if 860 end if 861 862 end subroutine setup_positron