TABLE OF CONTENTS


ABINIT/m_positron [ Modules ]

[ Top ] [ Modules ]

NAME

  m_positron

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_positron
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_efield
33  use m_errors
34  use m_abicore
35  use m_energies
36  use m_wffile
37  use m_electronpositron
38  use m_hdr
39  use m_xmpi
40  use m_bandfft_kpt
41 
42  use m_special_funcs,  only : sbf8
43  use m_ioarr,    only : ioarr, read_rhor
44  use m_pawang,   only : pawang_type, realgaunt
45  use m_pawrad,   only : pawrad_type, simp_gen, nderiv_gen
46  use m_pawtab,   only : pawtab_type
47  use m_paw_ij,   only : paw_ij_type
48  use m_pawfgrtab,only : pawfgrtab_type
49  use m_pawrhoij,only : pawrhoij_type, pawrhoij_copy, pawrhoij_alloc, pawrhoij_free,&
50                        pawrhoij_nullify, pawrhoij_gather, pawrhoij_get_nspden, symrhoij
51  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_send, &
52                         pawcprj_mpi_recv, pawcprj_free, pawcprj_copy, pawcprj_bcast
53  use m_pawfgr,   only : pawfgr_type
54  use m_paw_nhat, only : pawmknhat
55  use m_fock,     only : fock_type
56  use m_kg,       only : getcut
57  use defs_wvltypes,     only : wvl_data
58  use m_spacepar,        only : hartre
59  use m_mkrho,           only : initro
60  use m_paw_occupancies, only : initrhoij, pawaccrhoij
61  use m_gammapositron, only : gammapositron, gammapositron_fft
62  use m_forstr,          only : forstr
63  use m_pawxc,         only : pawxcsum
64  use m_paw_denpot,    only : pawdensities
65  use m_drivexc,       only : mkdenpos
66 
67  use m_paw_sphharm, only : initylmr
68  use m_pawpsp,  only : pawpsp_read_corewf
69  use m_crystal, only : crystal_t
70  use m_mpinfo,  only : ptabs_fourdp,set_mpi_enreg_fft,unset_mpi_enreg_fft,destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle
71  use m_io_tools,only : open_file,close_unit,get_unit
72  use m_fftcore, only : sphereboundary
73  use m_prep_kgb,        only : prep_fourwf
74  use m_fft,            only : fourwf, fourdp
75 
76  implicit none
77 
78  private

ABINIT/posdoppler [ Functions ]

[ Top ] [ 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
   | use_gpu_cuda=flag for Cuda use
   | 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= informations 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

PARENTS

      outscfcv

CHILDREN

      bandfft_kpt_destroy,bandfft_kpt_mpi_recv,bandfft_kpt_mpi_send
      bandfft_kpt_reset,destroy_mpi_enreg,fourdp,fourwf,gammapositron_fft
      initmpi_seq,initylmr,mkdenpos,pawaccrhoij,pawcprj_alloc,pawcprj_bcast
      pawcprj_copy,pawcprj_free,pawcprj_get,pawcprj_mpi_recv,pawcprj_mpi_send
      pawpsp_read_corewf,pawrhoij_alloc,pawrhoij_free,pawrhoij_gather
      pawrhoij_nullify,poslifetime,posratecore,prep_fourwf,ptabs_fourdp,sbf8
      set_mpi_enreg_fft,simp_gen,sphereboundary,symrhoij,unset_mpi_enreg_fft
      wffclose,wffopen,wrtout,xderivewrecend,xderivewrecinit,xderivewrite
      xmoveoff,xmpi_bcast,xmpi_recv,xmpi_send,xmpi_sum

SOURCE

1857 !Macro to go from row-column indexing to combined indexing
1858 #define RCC(glmn,hlmn) max(glmn,hlmn)*(max(glmn,hlmn)-1)/2+min(glmn,hlmn)
1859 !Macro to go from l,m angular momentum indexing to combined indexing
1860 #define LMC(lval,mval) lval*lval+lval+mval+1
1861 
1862 subroutine posdoppler(cg,cprj,Crystal,dimcprj,dtfil,dtset,electronpositron,&
1863 &                     filpsp,kg,mcg,mcprj,mpi_enreg,my_natom,&
1864 &                     n3xccc,nfft,ngfft,nhat,npwarr,occ,pawang,pawrad,&
1865 &                     pawrhoij,pawtab,rhor,xccc3d)
1866 
1867 
1868 !This section has been created automatically by the script Abilint (TD).
1869 !Do not modify the following lines by hand.
1870 #undef ABI_FUNC
1871 #define ABI_FUNC 'posdoppler'
1872 !End of the abilint section
1873 
1874  implicit none
1875 
1876 !Arguments ------------------------------------
1877 !scalars
1878  integer,intent(in) :: mcg,mcprj,my_natom,n3xccc,nfft
1879  type(crystal_t) :: Crystal
1880  type(datafiles_type),intent(in) :: dtfil
1881  type(dataset_type),intent(in) :: dtset
1882  type(electronpositron_type),pointer :: electronpositron
1883  type(MPI_type),intent(inout) :: mpi_enreg
1884  type(pawang_type),intent(in) :: pawang
1885 !arrays
1886  integer,intent(in) :: dimcprj(dtset%natom)
1887  integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),ngfft(18),npwarr(dtset%nkpt)
1888  real(dp),intent(in) :: nhat(nfft,dtset%nspden*dtset%usepaw),xccc3d(n3xccc)
1889  real(dp),intent(in),target :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
1890  real(dp),intent(in),target :: rhor(nfft,dtset%nspden)
1891  real(dp),intent(inout),target :: cg(2,mcg)
1892  character(len=fnlen),intent(in) :: filpsp(dtset%ntypat)
1893  type(pawcprj_type),target :: cprj(dtset%natom,mcprj)
1894  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw)
1895  type(pawrhoij_type),intent(in),target :: pawrhoij(my_natom*dtset%usepaw)
1896  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
1897 
1898 !Local variables-------------------------------
1899 !scalars
1900  integer :: accessfil,basis_size,bandpp,bdtot_index,bdtot_index_pos,blocksize,cplex,cplex_rhoij
1901  integer :: glmij,i0lmn,i1,i2,i3,iat,iatm,iatom
1902  integer :: ib,ib_cprj,ib_cprj_pos,ib_pos,ibg,ibg_pos
1903  integer :: iblock,iblock_pos,ibpp,ibpp_pos
1904  integer :: icg,icg_pos,id1,id2,id3,ierr,ig1,ig2,ig3,igamma,ii,ikg,ikg_pos,ikpt
1905  integer :: ikpt_pos,il,ilm,ilmn,iln,indx,indx0,iorder_cprj,iproc,ir,isppol,isppol_pos,istwf_k
1906  integer :: istwf_k_pos,itypat,iwarn,iwavef,iwavef_pos,j2,j3,jj,jkpt,jl,jlm,jlmn,jln
1907  integer :: klm,kln,klmn,l_size,ll,llmax,llmin,lm,lmn_size,lmn2_size
1908  integer :: mband_cprj,mband_cprj_pos,mcg_pos
1909  integer :: mcprj_k,mcprj_k_pos,me_band,me_fft,me_kpt,me_kptband
1910  integer :: mesh_size,meshsz,mm,my_ngrid,my_nspinor,my_nsppol,my_n2,n1,n2,n3,n4,n5,n6
1911  integer :: nband_cprj_eff_pos,nband_cprj_k,nband_cprj_k_pos
1912  integer :: nband_eff_pos,nband_k,nband_k_pos
1913  integer :: nblock_band,nblock_band_eff_pos,nkpt
1914  integer :: nproc_band,nproc_fft,nproc_kpt,nproc_kptband,npw_k,npw_k_pos
1915  integer :: nspden_rhoij,option,tag,unit_doppler
1916  integer :: tim_fourdp=0,tim_fourwf=-36
1917  integer :: ylmr_normchoice,ylmr_npts,ylmr_option
1918  logical,parameter :: include_nhat_in_gamma=.false.,state_dependent=.true.
1919  logical,parameter :: kgamma_only_positron=.true.,wf_conjugate=.false.
1920  logical :: cprj_paral_band,ex,mykpt,mykpt_pos,usetimerev
1921  real(dp) :: arg,bessarg,cpi,cpr,cp11,cp12,cp21,cp22,gammastate,intg
1922  real(dp) :: lambda_v1,lambda_v2,lambda_core,lambda_pw,occ_el,occ_pos
1923  real(dp) :: pnorm,pr,rate,rate_ipm,ratec,ratec_ipm,rate_paw,rate_paw_ipm
1924  real(dp) :: scale_,units_,weight,weight_pos,wf_fact,wtk_k,wtk_k_pos,vec
1925  character(len=fnlen) :: filename,filename_dop
1926  character(len=500) :: msg
1927  type(bandfft_kpt_type),pointer :: bandfft_kpt_el,bandfft_kpt_pos
1928  type(MPI_type) :: mpi_enreg_seq
1929  type(wffile_type) :: wff
1930 !arrays
1931  integer,allocatable :: gbound(:,:),gbound_pos(:,:),kg_k(:,:),kg_k_pos(:,:)
1932  integer,allocatable :: lcor(:),lmncmax(:),my_ffttab(:),my_gridtab(:),ncor(:),nphicor(:)
1933  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1934  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1935  logical,allocatable :: have_intc(:,:,:),have_rad(:,:)
1936  real(dp) :: buf(4),contrib(2),cp(2),cp_pos(2),expipr(2),pbn(3),pcart(3)
1937  real(dp) :: radsumnfftc(2),ylmgr(1,1,0),ylmr_nrm(1)
1938  real(dp),allocatable :: cwaveg(:,:),cwaveg_pos(:,:),cwaver(:),cwaver_pos(:),cwaver_pos_block(:)
1939  real(dp),allocatable :: cg_k_pos(:,:),cwaveaug(:,:,:,:),cwaveaug_pos(:,:,:,:)
1940  real(dp),allocatable :: denpot_dum(:,:,:),energycor(:),ff(:),fofgout_dum(:,:)
1941  real(dp),allocatable :: gamma(:,:),intc(:,:,:),j_bessel(:,:),jbes(:),mpibuf(:,:)
1942  real(dp),allocatable :: occ_k(:),occ_k_pos(:),pcart_k(:,:)
1943  real(dp),allocatable :: radint1(:,:),radint2(:,:),radint3(:,:)
1944  real(dp),allocatable :: radsumnfft1(:,:),radsumnfft2(:,:),radsumnfft3(:,:)
1945  real(dp),allocatable :: rho_contrib(:),rho_contrib_g(:,:)
1946  real(dp),allocatable :: rho_contrib_paw1(:,:),rho_contrib_paw2(:,:),rho_contrib_paw3(:,:)
1947  real(dp),allocatable :: rho_moment_v1(:,:),rho_moment_v2(:,:)
1948  real(dp),allocatable :: rho_moment_core(:,:),rho_moment_k(:),rho_moment_k2(:)
1949  real(dp),allocatable :: rho_pw(:,:),rhor_dop_el(:)
1950  real(dp),allocatable :: rhocorej(:),rhoe(:,:),rhop(:,:),ylmp(:)
1951  real(dp),pointer :: cg_pos_ptr(:,:),cg_ptr(:,:),occ_ptr(:),occ_pos_ptr(:)
1952  real(dp),pointer :: rhor_(:,:),rhor_ep_(:,:)
1953  complex(dpc) :: ifac ! (-i)^L mod 4
1954  complex(dpc),dimension(0:3) :: ilfac(0:3)=(/(1.0,0.0),(0.0,-1.0),(-1.0,0.0),(0.0,1.0)/)
1955  type(coeff1_type),allocatable :: gammastate_c(:)
1956  type(coeffi2_type),allocatable :: indlmncor(:)
1957  type(coeff2_type),allocatable :: phicor(:)
1958  type(coeff6_type),allocatable :: radsum1(:),radsum2(:),radsum3(:)
1959  type(coeff7_type),allocatable :: radsumc(:)
1960  type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_k_pos(:,:),cprj_pos(:,:)
1961  type(pawcprj_type),pointer :: cprj_pos_ptr(:,:),cprj_ptr(:,:)
1962  type(pawrhoij_type),allocatable :: pawrhoij_dop_el(:)
1963  type(pawrhoij_type),pointer :: pawrhoij_ptr(:),pawrhoij_all(:),pawrhoij_ep_all(:)
1964 
1965 ! *************************************************************************
1966 
1967  DBG_ENTER("COLL")
1968 
1969 !Compatibility tests
1970  if (.not.associated(electronpositron)) then
1971    msg='electronpositron variable must be associated!'
1972    MSG_BUG(msg)
1973  end if
1974  if (allocated(mpi_enreg%proc_distrb)) then
1975    do isppol=1,dtset%nsppol
1976      do ikpt=1,dtset%nkpt
1977        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1978        if (any(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)/=mpi_enreg%proc_distrb(ikpt,1,isppol))) then
1979          msg='proc_distrib cannot be distributed over bands!'
1980          MSG_BUG(msg)
1981        end if
1982      end do
1983    end do
1984  end if
1985  if (dtset%nspinor==2) then
1986    msg='Doppler broadening not available for spinorial wave functions (nspinor=2)!'
1987    MSG_BUG(msg)
1988  end if
1989  if (mcprj==0) then
1990    msg='<p|Psi> (cprj) datastructure must be kept in memory (see pawusecp input keyword)!'
1991    MSG_BUG(msg)
1992  end if
1993  if (dtset%usepaw==0) then
1994    write(msg,'(5a)') 'Momentum distribution of annihilating electron-positron pairs',ch10,&
1995 &   'in the Norm-conserving Pseudopotential formalism is incomplete!',ch10,&
1996 &   'No core contribution is included.'
1997    MSG_WARNING(msg)
1998  end if
1999  if (any(dtset%nband(:)/=dtset%nband(1))) then
2000    write(msg,'(a)') 'Number of bands has to be the same for all k-points!'
2001    MSG_BUG(msg)
2002  end if
2003  if (dtset%usepaw==1) then
2004    if (size(pawrhoij)/=mpi_enreg%my_natom) then
2005      write(msg,'(a)') 'wrong size for pawrhoij! '
2006      MSG_BUG(msg)
2007    end if
2008  end if
2009 
2010 !Various initializations
2011  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
2012  n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6)
2013  id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2
2014  iorder_cprj=0 ; cplex=2 ; iwarn=1
2015  wf_fact=one;if (wf_conjugate) wf_fact=-one
2016  nkpt=dtset%nkpt
2017 
2018 !Manage kpt/spin parallelism
2019  ABI_ALLOCATE(my_gridtab,(nkpt))
2020  my_gridtab=0
2021  do ii=1,nkpt
2022    if (any(mpi_enreg%my_isppoltab(:)==1)) my_gridtab(ii)=mpi_enreg%my_kpttab(ii)
2023  end do
2024  my_ngrid=count(my_gridtab(:)/=0)
2025  my_nsppol=sum(mpi_enreg%my_isppoltab(:))
2026  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
2027 
2028 !Parallel settings
2029  if (mpi_enreg%paral_kgb/=0) then
2030    nproc_kpt=mpi_enreg%nproc_kpt
2031    nproc_band=mpi_enreg%nproc_band
2032    nproc_fft=mpi_enreg%nproc_fft
2033    nproc_kptband=xmpi_comm_size(mpi_enreg%comm_kptband)
2034    me_kpt=mpi_enreg%me_kpt
2035    me_band=mpi_enreg%me_band
2036    me_fft=mpi_enreg%me_fft
2037    me_kptband=xmpi_comm_rank(mpi_enreg%comm_kptband)
2038    bandpp=mpi_enreg%bandpp
2039    my_n2=n2/nproc_fft
2040    accessfil=IO_MODE_FORTRAN;if(nproc_fft>1)accessfil=IO_MODE_MPI
2041  else
2042    nproc_kpt=mpi_enreg%nproc_kpt
2043    nproc_band=1;nproc_fft=1
2044    nproc_kptband=nproc_kpt
2045    me_band=0;me_fft=0
2046    me_kpt=mpi_enreg%me_kpt
2047    me_kptband=me_kpt
2048    bandpp=1 ; my_n2=n2
2049    accessfil=IO_MODE_FORTRAN
2050  end if
2051  blocksize=nproc_band*bandpp
2052  nblock_band=dtset%nband(1)/blocksize
2053 
2054 !Select density according to nhat choice0
2055  if (dtset%usepaw==0.or.include_nhat_in_gamma) then
2056    rhor_ => rhor
2057    rhor_ep_ => electronpositron%rhor_ep
2058  else
2059    ABI_ALLOCATE(rhor_,(nfft,dtset%nspden))
2060    ABI_ALLOCATE(rhor_ep_,(nfft,dtset%nspden))
2061    rhor_=rhor-nhat
2062    rhor_ep_=electronpositron%rhor_ep-electronpositron%nhat_ep
2063  end if
2064 
2065 !Select type(s) of enhancement factor
2066  igamma=0
2067  if (electronpositron%ixcpositron==-1) igamma=0
2068  if (electronpositron%ixcpositron== 1) igamma=2
2069  if (electronpositron%ixcpositron== 2) igamma=4
2070  if (electronpositron%ixcpositron== 3) igamma=2
2071  if (electronpositron%ixcpositron==11) igamma=3
2072  if (electronpositron%ixcpositron==31) igamma=3
2073 
2074 !Select electronic and positronic states
2075  if (electronpositron%particle==EP_ELECTRON) then !we should not be in this case
2076    cg_ptr => electronpositron%cg_ep
2077    cprj_ptr => electronpositron%cprj_ep
2078    occ_ptr => electronpositron%occ_ep
2079    pawrhoij_ptr => electronpositron%pawrhoij_ep
2080    cg_pos_ptr => cg
2081    cprj_pos_ptr => cprj
2082    occ_pos_ptr => occ
2083  end if
2084  if (electronpositron%particle==EP_POSITRON) then
2085    cg_ptr => cg
2086    cprj_ptr => cprj
2087    occ_ptr => occ
2088    pawrhoij_ptr => pawrhoij
2089    cg_pos_ptr => electronpositron%cg_ep
2090    cprj_pos_ptr => electronpositron%cprj_ep
2091    occ_pos_ptr => electronpositron%occ_ep
2092  end if
2093 
2094 !Determine if cprj datastructures are distributed over bands
2095  mband_cprj=size(cprj_ptr,2)/(my_nspinor*dtset%mkmem*dtset%nsppol)
2096  mband_cprj_pos=size(cprj_pos_ptr,2)/(my_nspinor*dtset%mkmem*dtset%nsppol)
2097  cprj_paral_band=(mband_cprj<dtset%mband)
2098 
2099 !Get the distrib associated with the fft_grid
2100  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2101 
2102 !===============================================================================
2103 !================ Calculate the PAW on-site constants ==========================
2104 
2105  if (dtset%usepaw==1) then
2106 
2107    ylmr_normchoice = 0 ! input to initylmr are normalized
2108    ylmr_npts = 1 ! only 1 point to compute in initylmr
2109    ylmr_nrm(1) = one ! weight of normed point for initylmr
2110    ylmr_option = 1 ! compute only ylm's in initylmr
2111 
2112   !Prepare radial integral for PAW correction for each atom type
2113    ABI_DATATYPE_ALLOCATE(radsum1,(dtset%ntypat))
2114    ABI_DATATYPE_ALLOCATE(radsum2,(dtset%ntypat))
2115    ABI_DATATYPE_ALLOCATE(radsum3,(dtset%ntypat))
2116    ABI_DATATYPE_ALLOCATE(radsumc,(dtset%ntypat))
2117 
2118    ABI_DATATYPE_ALLOCATE(indlmncor,(dtset%ntypat))
2119    ABI_DATATYPE_ALLOCATE(phicor,(dtset%ntypat))
2120    ABI_DATATYPE_ALLOCATE(gammastate_c,(dtset%natom))
2121    ABI_ALLOCATE(nphicor,(dtset%ntypat))
2122    ABI_ALLOCATE(lmncmax,(dtset%ntypat))
2123 
2124 !  Reading of core wave functions
2125    if (mpi_enreg%me_cell==0) then
2126      do itypat=1,dtset%ntypat
2127        filename=trim(filpsp(itypat))//'.corewf'
2128        inquire(file=filename,exist=ex)
2129        if (.not.ex) then
2130          write(unit=filename,fmt='(a,i1)') 'corewf.abinit',itypat
2131          inquire(file=filename,exist=ex)
2132          if (.not.ex) then
2133            msg='Core wave-functions file is missing!'
2134            MSG_ERROR(msg)
2135          end if
2136        end if
2137        call pawpsp_read_corewf(energycor,indlmncor(itypat)%value,lcor,lmncmax(itypat),&
2138 &       ncor,nphicor(itypat),pawrad(itypat),phicor(itypat)%value,&
2139 &       filename=filename)
2140 !      The following arrays are not used anymore
2141        ABI_DEALLOCATE(energycor)
2142        ABI_DEALLOCATE(lcor)
2143        ABI_DEALLOCATE(ncor)
2144      end do
2145    end if
2146    if (mpi_enreg%nproc_cell>1) then
2147      call xmpi_bcast(indlmncor,0,mpi_enreg%comm_cell,ierr)
2148      call xmpi_bcast(phicor,0,mpi_enreg%comm_cell,ierr)
2149      call xmpi_bcast(nphicor,0,mpi_enreg%comm_cell,ierr)
2150      call xmpi_bcast(lmncmax,0,mpi_enreg%comm_cell,ierr)
2151    end if
2152 
2153    do itypat=1,dtset%ntypat
2154 
2155      mesh_size = pawtab(itypat)%mesh_size
2156      l_size = pawtab(itypat)%l_size
2157      lmn_size = pawtab(itypat)%lmn_size
2158      lmn2_size = pawtab(itypat)%lmn2_size
2159      basis_size = pawtab(itypat)%basis_size
2160 
2161      ABI_ALLOCATE(j_bessel,(mesh_size,l_size))
2162      ABI_ALLOCATE(ylmp,(l_size*l_size))
2163      ABI_ALLOCATE(have_intc,(l_size,basis_size,nphicor(itypat)))
2164      ABI_ALLOCATE(have_rad,(l_size,pawtab(itypat)%ij_size))
2165      ABI_ALLOCATE(intc,(l_size,basis_size,nphicor(itypat)))
2166      ABI_ALLOCATE(radint1,(l_size,pawtab(itypat)%ij_size))
2167      ABI_ALLOCATE(radint2,(l_size,pawtab(itypat)%ij_size))
2168      ABI_ALLOCATE(radint3,(l_size,pawtab(itypat)%ij_size))
2169 
2170      ABI_ALLOCATE(radsumc(itypat)%value,(2,lmn_size,lmncmax(itypat),n1,my_n2,n3,my_ngrid))
2171      ABI_ALLOCATE(radsum1(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid))
2172      ABI_ALLOCATE(radsum2(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid))
2173      ABI_ALLOCATE(radsum3(itypat)%value,(2,lmn2_size,n1,my_n2,n3,my_ngrid))
2174      radsumc(itypat)%value=zero
2175      radsum1(itypat)%value=zero
2176      radsum2(itypat)%value=zero
2177      radsum3(itypat)%value=zero
2178 
2179      ABI_ALLOCATE(jbes,(l_size))
2180      ABI_ALLOCATE(ff,(mesh_size))
2181      meshsz=pawrad(itypat)%int_meshsz
2182      if (meshsz>mesh_size) ff(meshsz+1:mesh_size)=zero
2183 
2184      indx=0;jkpt=0
2185      do ikpt=1,nkpt
2186        if (my_gridtab(ikpt)==0) cycle
2187        jkpt=jkpt+1
2188        do i3=1,n3
2189          ig3=i3-(i3/id3)*n3-1
2190          do i2=1,n2
2191            if (me_fft/=fftn2_distrib(i2)) cycle
2192            j2=ffti2_local(i2)
2193            indx=n1*(my_n2*(i3-1)+(j2-1))
2194            ig2=i2-(i2/id2)*n2-1
2195            do i1=1,n1
2196              ig1=i1-(i1/id1)*n1-1
2197              indx=indx+1;if (mod(indx-1,nproc_band)/=me_band) cycle
2198 
2199              pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+&
2200 &             Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+&
2201 &             Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt))
2202              pnorm=dsqrt(dot_product(pcart,pcart))
2203 
2204              if (pnorm < tol12) then
2205                pbn(:) = zero
2206                ylmp(:) = zero
2207                ylmp(1) = 1.d0/sqrt(four_pi)
2208              else
2209                pbn(:) = pcart(:)/pnorm ! unit vector
2210                call initylmr(l_size,ylmr_normchoice,ylmr_npts,ylmr_nrm,ylmr_option,pbn(:),ylmp(:),ylmgr)
2211              end if
2212 
2213              pnorm=two_pi*pnorm ! re-normed for call to bessel
2214              do ir = 1, mesh_size
2215                bessarg = pnorm*pawrad(itypat)%rad(ir)
2216                call sbf8(l_size,bessarg,jbes)
2217                j_bessel(ir,:)=jbes(:)
2218              end do
2219 
2220 !            ===== Core part =====
2221 !            Need intc=\int phi phi_core jl (pr) dr
2222 
2223              have_intc(:,:,:)=.FALSE. ; intc(:,:,:)=zero
2224 
2225              do jlmn = 1,lmncmax(itypat)
2226                jln = indlmncor(itypat)%value(5,jlmn)
2227                jlm = indlmncor(itypat)%value(4,jlmn)
2228                jl  = indlmncor(itypat)%value(1,jlmn)
2229                do ilmn = 1,lmn_size
2230                  iln = pawtab(itypat)%indlmn(5,ilmn)
2231                  ilm = pawtab(itypat)%indlmn(4,ilmn)
2232                  il  = pawtab(itypat)%indlmn(1,ilmn)
2233 
2234                  llmin = abs(il-jl)
2235                  llmax = il+jl
2236                  klm = RCC(ilm,jlm)
2237                  do ll=llmin,llmax,2
2238                    ifac=ilfac(mod(ll,4))
2239 
2240                    if (.not.have_intc(ll+1,iln,jln)) then
2241                      ff(1:mesh_size)=(pawtab(itypat)%phi(1:mesh_size,iln)*phicor(itypat)%value(1:mesh_size,jln))&
2242 &                     *j_bessel(1:mesh_size,ll+1)
2243                      call simp_gen(intg,ff,pawrad(itypat))
2244                      intc(ll+1,iln,jln)=intg
2245                      have_intc(ll+1,iln,jln)=.true.
2246                    end if
2247 
2248                    do mm=-ll,ll
2249                      lm = LMC(ll,mm)
2250                      glmij=pawang%gntselect(lm,klm)
2251                      if (glmij>0) then
2252                        arg=ylmp(lm)*pawang%realgnt(glmij)*intc(ll+1,iln,jln)
2253                        radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt) = &
2254 &                       radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt)+arg*real(ifac)
2255                        radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt) = &
2256 &                       radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt)+arg*aimag(ifac)
2257                      end if
2258                    end do !mm
2259                  end do !ll
2260                end do !ilmn
2261              end do !jlmn
2262 
2263 !            ===== Valence part =====
2264 !            Need int1=\int phi_i phi_j jl (pr) dr
2265 !            and  int2=\int tphi_i tphi_j jl (pr) dr
2266 
2267              have_rad(:,:)= .FALSE.;radint1=zero;radint2=zero;radint3=zero
2268 
2269              do klmn=1,pawtab(itypat)%lmn2_size
2270                klm=pawtab(itypat)%indklmn(1,klmn);kln=pawtab(itypat)%indklmn(2,klmn)
2271                llmin=pawtab(itypat)%indklmn(3,klmn);llmax=pawtab(itypat)%indklmn(4,klmn)
2272 
2273                do ll=llmin,llmax,2
2274                  ifac=ilfac(mod(ll,4))
2275 
2276                  if (.not.have_rad(ll+1,kln)) then
2277                    ff(1:mesh_size)=pawtab(itypat)%phiphj(1:mesh_size,kln)*j_bessel(1:mesh_size,ll+1)
2278                    call simp_gen(intg,ff,pawrad(itypat))
2279                    radint1(ll+1,kln)=intg
2280                    ff(1:mesh_size)=pawtab(itypat)%tphitphj(1:mesh_size,kln)*j_bessel(1:mesh_size,ll+1)
2281                    call simp_gen(intg,ff,pawrad(itypat))
2282                    radint2(ll+1,kln)=intg
2283                    ff(1:mesh_size)=(pawtab(itypat)%phiphj  (1:mesh_size,kln) &
2284 &                   -pawtab(itypat)%tphitphj(1:mesh_size,kln))&
2285 &                   *j_bessel(1:mesh_size,ll+1)
2286                    call simp_gen(intg,ff,pawrad(itypat))
2287                    radint3(ll+1,kln)=intg
2288                    have_rad(ll+1,kln)=.true.
2289                  end if
2290 
2291                  do mm=-ll,ll
2292                    lm = LMC(ll,mm)
2293                    glmij=pawang%gntselect(lm,klm)
2294                    if (glmij>0) then
2295                      arg=ylmp(lm)*pawang%realgnt(glmij)
2296                      radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt) = &
2297 &                     radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint1(ll+1,kln)
2298                      radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt) = &
2299 &                     radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint1(ll+1,kln)
2300                      radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt) = &
2301 &                     radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint2(ll+1,kln)
2302                      radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt) = &
2303 &                     radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint2(ll+1,kln)
2304                      radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt) = &
2305 &                     radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt)+real(ifac) *arg*radint3(ll+1,kln)
2306                      radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt) = &
2307 &                     radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt)+aimag(ifac)*arg*radint3(ll+1,kln)
2308                    end if
2309                  end do !mm
2310                end do !ll
2311              end do !klmn
2312 
2313            end do ! end loop over i1
2314          end do ! end loop over i2
2315        end do ! end loop over i3
2316      end do ! end loop over ikpt
2317 
2318      ABI_DEALLOCATE(ff)
2319      ABI_DEALLOCATE(jbes)
2320 
2321      ABI_DEALLOCATE(j_bessel)
2322      ABI_DEALLOCATE(ylmp)
2323 
2324      ABI_DEALLOCATE(intc)
2325      ABI_DEALLOCATE(have_intc)
2326 
2327      ABI_DEALLOCATE(radint1)
2328      ABI_DEALLOCATE(radint2)
2329      ABI_DEALLOCATE(radint3)
2330      ABI_DEALLOCATE(have_rad)
2331 
2332      call xmpi_sum(radsumc(itypat)%value,mpi_enreg%comm_band,ierr)
2333      call xmpi_sum(radsum1(itypat)%value,mpi_enreg%comm_band,ierr)
2334      call xmpi_sum(radsum2(itypat)%value,mpi_enreg%comm_band,ierr)
2335      call xmpi_sum(radsum3(itypat)%value,mpi_enreg%comm_band,ierr)
2336 
2337    end do ! end loop over atom types
2338  end if ! PAW
2339 
2340 !Allocate main memory
2341  ABI_ALLOCATE(rho_contrib,(cplex*nfft))
2342  ABI_ALLOCATE(rho_contrib_g,(cplex,nfft))
2343  ABI_ALLOCATE(rho_contrib_paw1,(cplex,nfft))
2344  ABI_ALLOCATE(rho_contrib_paw2,(cplex,nfft))
2345  ABI_ALLOCATE(rho_contrib_paw3,(cplex,nfft))
2346 
2347  ABI_ALLOCATE(rho_moment_v1,(nfft,my_ngrid))
2348  ABI_ALLOCATE(rho_moment_v2,(nfft,my_ngrid))
2349  ABI_ALLOCATE(rho_moment_core,(nfft,my_ngrid))
2350  ABI_ALLOCATE(rho_pw,(nfft,my_ngrid))
2351  rho_moment_v1=zero;rho_moment_v2=zero
2352  rho_pw=zero;rho_moment_core=zero
2353 
2354 !Prepare gamma(r) for the state independent scheme
2355  ABI_ALLOCATE(gamma,(nfft,2))
2356  if (.not.state_dependent) then
2357    ABI_ALLOCATE(rhoe,(nfft,1))
2358    ABI_ALLOCATE(rhop,(nfft,1))
2359    if (electronpositron%particle==EP_ELECTRON) then
2360      rhoe(:,1)=rhor_ep_(:,1);rhop(:,1)=rhor_(:,1)
2361    else if (electronpositron%particle==EP_POSITRON) then
2362      rhoe(:,1)=rhor_(:,1);rhop(:,1)=rhor_ep_(:,1)
2363    end if
2364    call mkdenpos(iwarn,nfft,1,1,rhoe(:,1),dtset%xc_denpos)
2365    call mkdenpos(iwarn,nfft,1,1,rhop(:,1),dtset%xc_denpos)
2366    call gammapositron_fft(electronpositron,gamma,Crystal%gprimd,igamma,mpi_enreg,&
2367 &   n3xccc,nfft,ngfft,rhoe(:,1),rhop(:,1),xccc3d)
2368    ABI_DEALLOCATE(rhoe)
2369    ABI_DEALLOCATE(rhop)
2370  else
2371    gamma=one
2372  end if
2373 
2374 !Some allocations for state-dependent scheme
2375  if (state_dependent) then
2376 !  Fake MPI data to be used in poslifetime; allow only FFT parallelism
2377    call initmpi_seq(mpi_enreg_seq)
2378    mpi_enreg_seq%my_natom=dtset%natom
2379    call set_mpi_enreg_fft(mpi_enreg_seq,mpi_enreg%comm_fft,mpi_enreg%distribfft,&
2380 &   mpi_enreg%me_g0,mpi_enreg%paral_kgb)
2381 !  Allocate memory for state-dependent scheme
2382    ABI_ALLOCATE(rhor_dop_el,(nfft))
2383    if (dtset%usepaw==1) then
2384      ABI_DATATYPE_ALLOCATE(pawrhoij_dop_el,(dtset%natom))
2385      nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
2386      call pawrhoij_alloc(pawrhoij_dop_el,dtset%pawcpxocc,nspden_rhoij,&
2387      dtset%nspinor,dtset%nsppol,dtset%typat,&
2388      pawtab=pawtab,use_rhoij_=1,use_rhoijp=1)
2389 !    Cancel distribution of PAW data over atomic sites
2390 !    We use here pawrhoij because polifetime routine
2391 !    detects by itself the particle described by pawrhoij
2392      if (mpi_enreg%my_natom<dtset%natom) then
2393        ABI_DATATYPE_ALLOCATE(pawrhoij_all,(dtset%natom))
2394        call pawrhoij_nullify(pawrhoij_all)
2395        call pawrhoij_gather(pawrhoij,pawrhoij_all,-1,mpi_enreg%comm_atom, &
2396 &       with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.)
2397        ABI_DATATYPE_ALLOCATE(pawrhoij_ep_all,(dtset%natom))
2398        call pawrhoij_nullify(pawrhoij_ep_all)
2399        call pawrhoij_gather(electronpositron%pawrhoij_ep,pawrhoij_ep_all,-1,mpi_enreg%comm_atom, &
2400 &       with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.)
2401      else
2402        pawrhoij_all => pawrhoij
2403        pawrhoij_ep_all => electronpositron%pawrhoij_ep
2404      end if
2405    end if
2406  end if
2407 
2408 !==============================================================================
2409 !================ Loop over positronic states =================================
2410 
2411 !LOOP OVER k POINTS
2412  ibg_pos=0;icg_pos=0;ikg_pos=0;bdtot_index_pos=0;isppol_pos=1
2413  do ikpt_pos=1,merge(1,nkpt,kgamma_only_positron)
2414 
2415 !  Extract data for this kpt_pos
2416    npw_k_pos=npwarr(ikpt_pos)
2417    wtk_k_pos=dtset%wtk(ikpt_pos); if (kgamma_only_positron) wtk_k_pos=one
2418    istwf_k_pos=dtset%istwfk(ikpt_pos)
2419    nband_k_pos=dtset%nband(ikpt_pos+(isppol_pos-1)*nkpt)
2420    nband_cprj_k_pos=nband_k_pos/nproc_band
2421    mykpt_pos=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_pos,1,nband_k_pos,&
2422 &   isppol_pos,mpi_enreg%me_kpt))
2423 
2424 !  Retrieve additional data for this kpt_pos
2425    ABI_ALLOCATE(occ_k_pos,(nband_k_pos))
2426    occ_k_pos(:)=occ_pos_ptr(1+bdtot_index_pos:nband_k_pos+bdtot_index_pos)
2427    nband_eff_pos=1
2428    do ib_pos=1,nband_k_pos
2429      if (occ_k_pos(ib_pos)>tol8) nband_eff_pos=ib_pos
2430    end do
2431    if (mod(nband_eff_pos,blocksize)/=0) nband_eff_pos=((nband_eff_pos/blocksize)+1)*blocksize
2432 
2433    nblock_band_eff_pos=nband_eff_pos/blocksize
2434 
2435    mcg_pos=npw_k_pos*my_nspinor*nband_eff_pos
2436    ABI_ALLOCATE(cg_k_pos,(2,mcg_pos))
2437 
2438    mcprj_k_pos=0
2439    if (dtset%usepaw==1) then
2440      nband_cprj_eff_pos=nband_eff_pos/nproc_band
2441      mcprj_k_pos=my_nspinor*nband_cprj_eff_pos
2442      ABI_DATATYPE_ALLOCATE(cprj_k_pos,(dtset%natom,mcprj_k_pos))
2443      call pawcprj_alloc(cprj_k_pos,0,dimcprj)
2444    end if
2445 
2446    if (mpi_enreg%paral_kgb==0) then
2447      ABI_ALLOCATE(gbound_pos,(2*dtset%mgfft+8,2))
2448      ABI_ALLOCATE(kg_k_pos,(3,npw_k_pos))
2449    else if (mykpt_pos) then
2450      nullify(bandfft_kpt_pos)
2451    else
2452      ABI_DATATYPE_ALLOCATE(bandfft_kpt_pos,)
2453      call bandfft_kpt_reset(bandfft_kpt_pos)
2454    end if
2455 
2456 !  Exchange data (WF components) between procs
2457    if (mykpt_pos) then
2458      cg_k_pos(:,1:mcg_pos)=cg_pos_ptr(:,icg_pos+1:icg_pos+mcg_pos)
2459      if (mpi_enreg%paral_kgb==0) kg_k_pos(:,1:npw_k_pos)=kg(:,1+ikg_pos:npw_k_pos+ikg_pos)
2460      if (dtset%usepaw==1) then
2461        call pawcprj_get(Crystal%atindx1,cprj_k_pos,cprj_pos_ptr,dtset%natom,1,ibg_pos,ikpt_pos,iorder_cprj,&
2462 &       isppol_pos,mband_cprj_pos,dtset%mkmem,dtset%natom,nband_cprj_eff_pos,nband_k_pos,my_nspinor,&
2463 &       dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2464      end if
2465      if (mpi_enreg%paral_kgb/=0) then
2466        jj=mpi_enreg%my_kpttab(ikpt_pos)
2467        bandfft_kpt_pos => bandfft_kpt(jj)
2468      end if
2469      do ii=0,mpi_enreg%nproc_kpt-1
2470        if (ii/=mpi_enreg%me_kpt) then
2471          tag=ikpt_pos+(isppol_pos-1)*nkpt+2*nkpt*ii
2472          call xmpi_send(cg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2473          tag=tag+nkpt*(1+2*mpi_enreg%nproc_kpt)
2474          if (mpi_enreg%paral_kgb==0) then
2475            call xmpi_send(kg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2476          else
2477            call bandfft_kpt_mpi_send(bandfft_kpt_pos,ii,tag,mpi_enreg%comm_kpt,ierr,profile='fourwf')
2478          end if
2479          if (dtset%usepaw==1) then
2480            call pawcprj_mpi_send(dtset%natom,mcprj_k_pos,dimcprj,0,cprj_k_pos,ii,mpi_enreg%comm_kpt,ierr)
2481          end if
2482        end if
2483      end do
2484    else
2485      ii=0;if (allocated(mpi_enreg%proc_distrb)) ii=mpi_enreg%proc_distrb(ikpt_pos,1,isppol_pos)
2486      tag=ikpt_pos+(isppol_pos-1)*nkpt+2*nkpt*mpi_enreg%me_kpt
2487      call xmpi_recv(cg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2488      tag=tag+nkpt*(1+2*mpi_enreg%nproc_kpt)
2489      if (mpi_enreg%paral_kgb==0) then
2490        call xmpi_recv(kg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2491      else
2492        call bandfft_kpt_mpi_recv(bandfft_kpt_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
2493      end if
2494      if (dtset%usepaw==1) then
2495        call pawcprj_mpi_recv(dtset%natom,mcprj_k_pos,dimcprj,0,cprj_k_pos,ii,mpi_enreg%comm_kpt,ierr)
2496      end if
2497    end if
2498 
2499    if (mpi_enreg%paral_kgb==0) then
2500      call sphereboundary(gbound_pos,istwf_k_pos,kg_k_pos,dtset%mgfft,npw_k_pos)
2501    end if
2502 
2503    ABI_ALLOCATE(cwaver_pos,(cplex*nfft))
2504    ABI_ALLOCATE(cwaver_pos_block,(cplex*nfft*bandpp))
2505    if (dtset%usepaw==1) then
2506      ABI_DATATYPE_ALLOCATE(cprj_pos,(dtset%natom,my_nspinor))
2507      call pawcprj_alloc(cprj_pos,0,dimcprj)
2508    end if
2509 
2510 !  ============================================================================
2511 !  Loops on positronic bands
2512 
2513    do iblock_pos=1,nblock_band_eff_pos
2514      ib_pos=1+(iblock_pos-1)*blocksize
2515      if (any(abs(occ_k_pos(ib_pos:ib_pos+blocksize-1))>tol8)) then
2516 
2517        ABI_ALLOCATE(cwaveg_pos,(2,npw_k_pos*blocksize))
2518        ABI_ALLOCATE(cwaveaug_pos,(2,n4,n5,n6*bandpp))
2519        ABI_ALLOCATE(denpot_dum,(n4,n5,n6))
2520        ABI_ALLOCATE(fofgout_dum,(2,npw_k_pos*blocksize))
2521        iwavef_pos=(iblock_pos-1)*npw_k_pos*blocksize
2522        cwaveg_pos(:,1:npw_k_pos*blocksize)= &
2523 &       cg_k_pos(:,iwavef_pos+1:iwavef_pos+npw_k_pos*blocksize)
2524 
2525 !      Get positronic wave function in real space
2526        option=0
2527        if (mpi_enreg%paral_kgb==0) then
2528          weight_pos=occ_k_pos(ib_pos)*wtk_k_pos
2529          call fourwf(1,denpot_dum,cwaveg_pos,fofgout_dum,cwaveaug_pos,&
2530 &         gbound_pos,gbound_pos,istwf_k_pos,kg_k_pos,kg_k_pos,&
2531 &         dtset%mgfft,mpi_enreg,1,ngfft,npw_k_pos,npw_k_pos,&
2532 &         n4,n5,n6,option,mpi_enreg%paral_kgb,tim_fourwf,weight_pos,weight_pos,&
2533 &         use_gpu_cuda=dtset%use_gpu_cuda)
2534        else
2535          call prep_fourwf(denpot_dum,blocksize,cwaveg_pos,cwaveaug_pos,&
2536 &         iblock_pos,istwf_k_pos,dtset%mgfft,mpi_enreg,nband_k_pos,&
2537 &         bandpp,ngfft,npw_k_pos,n4,n5,n6,occ_k_pos,option,Crystal%ucvol,wtk_k_pos,&
2538 &         bandfft_kpt_tab=bandfft_kpt_pos,use_gpu_cuda=dtset%use_gpu_cuda)
2539        end if
2540 
2541        cwaver_pos_block=zero
2542        do ii=1,bandpp
2543          j3=(ii-1)*n3
2544          indx0=1+(ii-1)*cplex*nfft
2545          do i3=1,n3
2546            if (me_fft==fftn3_distrib(i3)) then
2547              indx=indx0+cplex*n1*n2*(ffti3_local(i3)-1)
2548              do i2=1,n2
2549                do i1=1,n1
2550                  cwaver_pos_block(indx  )=cwaveaug_pos(1,i1,i2,i3+j3)
2551                  cwaver_pos_block(indx+1)=cwaveaug_pos(2,i1,i2,i3+j3)
2552                  indx=indx+2
2553                end do
2554              end do
2555            end if
2556          end do
2557        end do
2558        ABI_DEALLOCATE(fofgout_dum)
2559        ABI_DEALLOCATE(denpot_dum)
2560        ABI_DEALLOCATE(cwaveaug_pos)
2561        ABI_DEALLOCATE(cwaveg_pos)
2562 
2563 !      At this stage, each band proc has bandpp bands in real space
2564 !      (distributed on FFT procs)
2565 
2566 !      ========================================================================
2567 !      Compute core contribution for this positronic band (PAW only)
2568 
2569        if (dtset%usepaw==1) then
2570          do ibpp_pos=1,bandpp
2571            ib_cprj_pos=(iblock_pos-1)*bandpp+ibpp_pos
2572            weight_pos=occ_k_pos(ib_pos+ibpp_pos-1+me_band*bandpp)*wtk_k_pos
2573 !       Calculate the annihilation rate for each core state for state dependent scheme
2574            iatm=0
2575            do itypat=1,dtset%ntypat
2576              mesh_size = pawtab(itypat)%mesh_size
2577              do iat=1,Crystal%nattyp(itypat)
2578                iatm=iatm+1;iatom=Crystal%atindx1(iatm)
2579                ABI_ALLOCATE(gammastate_c(iatom)%value,(lmncmax(itypat)))
2580                do jlmn=1,lmncmax(itypat)
2581                  jln = indlmncor(itypat)%value(5,jlmn)
2582                  contrib(:)=zero
2583                  ABI_ALLOCATE(rhocorej,(mesh_size))
2584                  rhocorej(1:mesh_size)=2*phicor(itypat)%value(1:mesh_size,jln)**2
2585                  call posratecore(dtset,electronpositron,iatom,dtset%natom,mesh_size,mpi_enreg_seq,&
2586 &                 1,pawang,pawrad,pawrhoij_all,pawrhoij_ep_all,pawtab,ratec,rhocorej)
2587 
2588                  call posratecore(dtset,electronpositron,iatom,dtset%natom,mesh_size,mpi_enreg_seq,&
2589 &                 2,pawang,pawrad,pawrhoij_all,pawrhoij_ep_all,pawtab,ratec_ipm,rhocorej)
2590 
2591                  gammastate_c(iatom)%value(jlmn)=ratec/ratec_ipm
2592                  ABI_DEALLOCATE(rhocorej)
2593                end do
2594              end do
2595            end do
2596            jkpt=0
2597            do ikpt=1,nkpt
2598              if (my_gridtab(ikpt)==0) cycle
2599              jkpt=jkpt+1
2600              do i3=1,n3
2601                ig3=i3-(i3/id3)*n3-1
2602                do i2=1,n2
2603                  if (me_fft==fftn2_distrib(i2)) then
2604                    j2=ffti2_local(i2)
2605                    ig2=i2-(i2/id2)*n2-1
2606                    indx=n1*(my_n2*(i3-1)+(j2-1))
2607                    do i1=1,n1
2608                      ig1=i1-(i1/id1)*n1-1
2609                      indx=indx+1
2610 
2611 !                    Loop on atoms (type sorted)
2612                      iatm=0
2613                      do itypat=1,dtset%ntypat
2614                        lmn_size = pawtab(itypat)%lmn_size
2615 
2616                        do iat=1,Crystal%nattyp(itypat)
2617                          iatm=iatm+1;iatom=Crystal%atindx1(iatm)
2618 
2619                          pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+&
2620 &                         Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+&
2621 &                         Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt))
2622                          pnorm=dsqrt(dot_product(pcart,pcart))
2623                          pr=dot_product(pcart,Crystal%xcart(:,iatom))
2624                          expipr(1)= cos(two_pi*pr)
2625                          expipr(2)=-sin(two_pi*pr)
2626 
2627 !                        Loop on ij states
2628                          do jlmn = 1,lmncmax(itypat)
2629                            contrib(:)=zero
2630                            do ilmn = 1,lmn_size
2631                              radsumnfftc(1)=expipr(1)*radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt)&
2632 &                             -expipr(2)*radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt)
2633                              radsumnfftc(2)=expipr(1)*radsumc(itypat)%value(2,ilmn,jlmn,i1,j2,i3,jkpt)&
2634 &                             +expipr(2)*radsumc(itypat)%value(1,ilmn,jlmn,i1,j2,i3,jkpt)
2635                              cp_pos(:)=cprj_k_pos(iatom,ib_cprj_pos)%cp(:,ilmn)
2636                              contrib(1)=contrib(1)+four_pi*(cp_pos(1)*radsumnfftc(1) &
2637 &                             -cp_pos(2)*radsumnfftc(2))
2638                              contrib(2)=contrib(2)+four_pi*(cp_pos(1)*radsumnfftc(2) &
2639 &                             +cp_pos(2)*radsumnfftc(1))
2640                            end do ! end loop over ilmn
2641                            ! 2 - electron state weight for 2 spins
2642                            rho_moment_core(indx,jkpt) = rho_moment_core(indx,jkpt) &
2643 &                           +gammastate_c(iatom)%value(jlmn)*2*weight_pos*(contrib(1)**2+contrib(2)**2)
2644                          end do ! end loop over jlmn
2645 
2646                        end do !end loop over atoms
2647                      end do !end loop over atom types
2648 
2649                    end do ! end loop over i1
2650                  end if ! end loop over i2
2651                end do
2652              end do ! end loop over i3
2653            end do ! jkpt
2654          end do ! ibpp_pos
2655        end if
2656 
2657 !      We now loop over positronic bands inside a block
2658 !      and select occupied ones
2659        do ibpp_pos=1,blocksize
2660          ib_pos=(iblock_pos-1)*blocksize+ibpp_pos
2661          occ_pos=occ_k_pos(ib_pos)
2662          if (abs(occ_pos)>tol8) then
2663 
2664 !          Parallelism: dirty trick (broadcast bands) but there should be few positronic bands (~1)
2665            if (nproc_band>1) then
2666              iproc=(ibpp_pos-1)/bandpp
2667              if (me_band==iproc) then
2668                indx=mod((ibpp_pos-1),bandpp)*cplex*nfft
2669                cwaver_pos(1:cplex*nfft)=cwaver_pos_block(indx+1:indx+cplex*nfft)
2670              end if
2671              call xmpi_bcast(cwaver_pos,iproc,mpi_enreg%comm_band,ierr)
2672              if (dtset%usepaw==1) then
2673                if (me_band==iproc) then
2674                  indx=mod((ibpp_pos-1),bandpp)*my_nspinor
2675                  call pawcprj_copy(cprj_k_pos(:,indx+1:indx+my_nspinor),cprj_pos)
2676                end if
2677                call pawcprj_bcast(cprj_pos,dtset%natom,my_nspinor,dimcprj,0,iproc,&
2678 &               mpi_enreg%comm_band,ierr)
2679              end if
2680            else
2681              cwaver_pos(1:cplex*nfft)=cwaver_pos_block(1:cplex*nfft)
2682              if (dtset%usepaw==1) then
2683                call pawcprj_copy(cprj_k_pos(:,(ib_pos-1)*my_nspinor+1:ib_pos*my_nspinor),cprj_pos)
2684              end if
2685            end if
2686 
2687 !      ========================================================================
2688 !      ================ Loop over electronic states ===========================
2689 
2690 !          Loop over spins
2691            ibg=0;icg=0;ikg=0;bdtot_index=0
2692            do isppol=1,dtset%nsppol
2693 !            Loop over k points
2694              ikg=0;jkpt=0
2695              do ikpt=1,nkpt
2696 
2697 !              Extract data for this kpt_pos
2698                npw_k=npwarr(ikpt)
2699                wtk_k=dtset%wtk(ikpt)
2700                istwf_k=dtset%istwfk(ikpt)
2701                nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
2702                nband_cprj_k=nband_k/nproc_band
2703                mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,&
2704 &               isppol,mpi_enreg%me_kpt))
2705 
2706 !              Select k-points for current proc
2707                if (mykpt) then
2708 
2709 !                Retrieve additional data for this kpt_pos
2710                  jkpt=jkpt+1
2711                  ABI_ALLOCATE(occ_k,(nband_k))
2712                  occ_k(:)=occ_ptr(1+bdtot_index:nband_k+bdtot_index)
2713 
2714                  mcprj_k=0
2715                  if (dtset%usepaw==1) then
2716                    mcprj_k=my_nspinor*nband_cprj_k
2717                    ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natom,mcprj_k))
2718                    call pawcprj_alloc(cprj_k,0,dimcprj)
2719                    call pawcprj_get(Crystal%atindx1,cprj_k,cprj_ptr,dtset%natom,1,ibg,ikpt,iorder_cprj,&
2720 &                   isppol,mband_cprj,dtset%mkmem,dtset%natom,nband_cprj_k,nband_cprj_k,my_nspinor,&
2721 &                   dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2722                  end if
2723 
2724                  if (mpi_enreg%paral_kgb==0) then
2725                    ABI_ALLOCATE(gbound,(2*dtset%mgfft+8,2))
2726                    ABI_ALLOCATE(kg_k,(3,npw_k))
2727                    kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
2728                    call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k)
2729                  else
2730                    jj=mpi_enreg%my_kpttab(ikpt)
2731                    bandfft_kpt_el => bandfft_kpt(jj)
2732                  end if
2733 
2734                  ABI_ALLOCATE(cwaver,(cplex*nfft*bandpp))
2735 
2736 !                ==================================================================
2737 !                Loops on electronic bands
2738 
2739                  do iblock=1,nblock_band
2740                    ib=1+(iblock-1)*blocksize
2741 
2742                    if (any(abs(occ_k(ib:ib+blocksize-1))>tol8)) then
2743 
2744 !                    Retrieve electronic wave function
2745                      ABI_ALLOCATE(cwaveg,(2,npw_k*blocksize))
2746                      ABI_ALLOCATE(cwaveaug,(2,n4,n5,n6*bandpp))
2747                      ABI_ALLOCATE(denpot_dum,(n4,n5,n6))
2748                      ABI_ALLOCATE(fofgout_dum,(2,npw_k*blocksize))
2749                      iwavef=(iblock-1)*npw_k*blocksize
2750                      cwaveg(:,1:npw_k*blocksize)= &
2751 &                     cg_ptr(:,icg+iwavef+1:icg+iwavef+npw_k*blocksize)
2752 
2753 !                    Get electronic wave function in real space
2754                      option=0
2755                      if (mpi_enreg%paral_kgb==0) then
2756                        weight=occ_k(ib)*wtk_k
2757                        call fourwf(1,denpot_dum,cwaveg,fofgout_dum,cwaveaug,&
2758 &                       gbound,gbound,istwf_k,kg_k,kg_k,&
2759 &                       dtset%mgfft,mpi_enreg,1,ngfft,npw_k,npw_k,&
2760 &                       n4,n5,n6,option,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,&
2761 &                       use_gpu_cuda=dtset%use_gpu_cuda)
2762                      else
2763                        call prep_fourwf(denpot_dum,blocksize,cwaveg,cwaveaug,&
2764 &                       iblock,istwf_k,dtset%mgfft,mpi_enreg,nband_k,&
2765 &                       bandpp,ngfft,npw_k,n4,n5,n6,occ_k,option,Crystal%ucvol,wtk_k,&
2766 &                       bandfft_kpt_tab=bandfft_kpt_el,use_gpu_cuda=dtset%use_gpu_cuda)
2767                      end if
2768 
2769                      cwaver=zero
2770                      do ii=1,bandpp
2771                        j3=(ii-1)*n3
2772                        indx0=1+(ii-1)*cplex*nfft
2773                        do i3=1,n3
2774                          if (me_fft==fftn3_distrib(i3)) then
2775                            indx=indx0+cplex*n1*n2*(ffti3_local(i3)-1)
2776                            do i2=1,n2
2777                              do i1=1,n1
2778                                cwaver(indx  )=cwaveaug(1,i1,i2,i3+j3)
2779                                cwaver(indx+1)=cwaveaug(2,i1,i2,i3+j3)
2780                                indx=indx+2
2781                              end do
2782                            end do
2783                          end if
2784                        end do
2785                      end do
2786                      ABI_DEALLOCATE(fofgout_dum)
2787                      ABI_DEALLOCATE(denpot_dum)
2788                      ABI_DEALLOCATE(cwaveaug)
2789                      ABI_DEALLOCATE(cwaveg)
2790 !                    At this stage, each band proc has bandpp bands in real space
2791 !                   (distributed on FFT procs)
2792 
2793 !                    We now loop on the bandpp bands
2794 !                    and select occupied ones
2795                      do ibpp=1,bandpp
2796                        occ_el=occ_k(ib+ibpp-1+me_band*bandpp)
2797                        if (abs(occ_el)>tol8) then
2798 
2799 !                        ==============================================================
2800 !                        Compute state-dependent annihilation rate
2801 !                        Avoid parallelism over kpt/bands/atoms
2802                          gammastate=one;rate_paw=one
2803                          if (state_dependent) then
2804                            weight=occ_el*wtk_k
2805                            ib_cprj=(iblock-1)*bandpp+ibpp
2806                            indx=1+(ibpp-1)*cplex*nfft
2807                            do ii=1,nfft
2808                              rhor_dop_el(ii)=weight*(cwaver(indx)*cwaver(indx)+cwaver(indx+1)*cwaver(indx+1))
2809                              indx=indx+2
2810                            end do
2811                            if (dtset%usepaw==1) then
2812                              do iatom=1,dtset%natom
2813                                pawrhoij_dop_el(iatom)%rhoij_=zero
2814                              end do
2815                              cplex_rhoij=2;if (istwf_k>1) cplex_rhoij=1
2816                              usetimerev=(dtset%kptopt>0.and.dtset%kptopt<3)
2817                              call pawaccrhoij(Crystal%atindx,cplex_rhoij,cprj_k(:,ib_cprj),cprj_k(:,ib_cprj),0,isppol,&
2818 &                             dtset%natom,dtset%natom,dtset%nspinor,occ_el,1,pawrhoij_dop_el,usetimerev,wtk_k)
2819 !                            Is it correct to apply symetries here (on a single band)?
2820 !                            If not, call symrhoij with nsym=1
2821                              call symrhoij(pawrhoij_dop_el,pawrhoij_dop_el,1,Crystal%gprimd,Crystal%indsym,0,dtset%natom,&
2822 &                             Crystal%nsym,dtset%ntypat,1,pawang,-10001,pawtab,Crystal%rprimd,Crystal%symafm,&
2823 &                             Crystal%symrec,dtset%typat)
2824                            end if
2825 !                          Has to call poslifetime in sequential because we are in a parallel section
2826 !                          Only FFT parallelism is allowed
2827                            call poslifetime(dtset,electronpositron,Crystal%gprimd,dtset%natom,mpi_enreg_seq,n3xccc,&
2828 &                           nfft,ngfft,nhat,2,pawang,pawrad,pawrhoij_all,pawtab,rate,rate_paw,rhor,Crystal%ucvol,xccc3d,&
2829 &                           rhor_dop_el=rhor_dop_el,pawrhoij_dop_el=pawrhoij_dop_el,pawrhoij_ep=pawrhoij_ep_all)
2830                            call poslifetime(dtset,electronpositron,Crystal%gprimd,dtset%natom,mpi_enreg_seq,n3xccc,&
2831 &                           nfft,ngfft,nhat,3,pawang,pawrad,pawrhoij_all,pawtab,rate_ipm,rate_paw_ipm,rhor,Crystal%ucvol,xccc3d,&
2832 &                           rhor_dop_el=rhor_dop_el,pawrhoij_dop_el=pawrhoij_dop_el,pawrhoij_ep=pawrhoij_ep_all)
2833                            gammastate=rate/rate_ipm
2834                            rate_paw=rate_paw/rate_paw_ipm
2835                          end if
2836 
2837 !                        ==============================================================
2838 !                        Compute plane-wave contribution to momentum distribution
2839 
2840 !                        Compute Psi^+(r) * Psi^-(r) * gamma(r) in real space
2841                          rho_contrib(:)=zero
2842                          indx=(ibpp-1)*cplex*nfft
2843                          if (cplex==2) then
2844                            do jj=1,nfft
2845                              ii=2*jj-1
2846                              rho_contrib(ii)  =sqrt(gamma(jj,2))*(cwaver_pos(ii)*cwaver(indx+ii)&
2847 &                             -wf_fact*cwaver_pos(ii+1)*cwaver(indx+ii+1))
2848                              rho_contrib(ii+1)=sqrt(gamma(jj,2))*(cwaver_pos(ii)*cwaver(indx+ii+1) &
2849 &                             +wf_fact*cwaver_pos(ii+1)*cwaver(indx+ii))
2850                            end do
2851                          else
2852                            do ii=1,nfft
2853                              rho_contrib(ii)=sqrt(gamma(ii,2))*cwaver_pos(ii)*cwaver(indx+ii)
2854                            end do
2855                          end if
2856 
2857 !                        FFT of (Psi+.Psi-.gamma) to get Intg[(Psi+.Psi-.gamma).exp(-igr)]
2858                          call fourdp(cplex,rho_contrib_g,rho_contrib,-1,mpi_enreg,nfft,ngfft,&
2859 &                         mpi_enreg%paral_kgb,tim_fourdp)
2860 
2861                          rho_pw(1:nfft,jkpt)=rho_pw(1:nfft,jkpt) +gammastate*occ_el*occ_pos &
2862 &                         *(rho_contrib_g(1,1:nfft)**2+rho_contrib_g(2,1:nfft)**2)
2863 
2864 !                        ==============================================================
2865 !                        Compute PAW on-site contribution to momentum distribution
2866 
2867                          if (dtset%usepaw==1) then
2868 
2869                            rho_contrib_paw1(:,:)= zero
2870                            rho_contrib_paw2(:,:)= zero
2871                            rho_contrib_paw3(:,:)= zero
2872 
2873                            ib_cprj=(iblock-1)*bandpp+ibpp
2874 
2875 !                          Loop on moments
2876                            indx=0
2877                            do i3=1,n3
2878                              ig3=i3-(i3/id3)*n3-1
2879                              do i2=1,n2
2880                                if (me_fft==fftn2_distrib(i2)) then
2881                                  j2=ffti2_local(i2)
2882                                  ig2=i2-(i2/id2)*n2-1
2883                                  indx=n1*(my_n2*(i3-1)+(j2-1))
2884                                  do i1=1,n1
2885                                    ig1=i1-(i1/id1)*n1-1
2886                                    indx=indx+1
2887 
2888                                    pcart(:)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+&
2889 &                                   Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+&
2890 &                                   Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt))
2891                                    pnorm=dsqrt(dot_product(pcart,pcart))
2892 
2893 !                                  Loop on atoms (type-sorted)
2894                                    iatm=0
2895                                    do itypat=1,dtset%ntypat
2896                                      lmn_size=pawtab(itypat)%lmn_size
2897                                      lmn2_size=pawtab(itypat)%lmn2_size
2898                                      ABI_ALLOCATE(radsumnfft1,(2,lmn2_size))
2899                                      ABI_ALLOCATE(radsumnfft2,(2,lmn2_size))
2900                                      ABI_ALLOCATE(radsumnfft3,(2,lmn2_size))
2901 
2902                                      do iat=1,Crystal%nattyp(itypat)
2903                                        iatm=iatm+1;iatom=Crystal%atindx1(iatm)
2904 
2905                                        pr=dot_product(pcart,Crystal%xcart(:,iatom))
2906                                        expipr(1)= cos(two_pi*pr)
2907                                        expipr(2)=-sin(two_pi*pr)
2908 
2909                                        do klmn=1,lmn2_size
2910                                          radsumnfft1(1,klmn)=expipr(1)*radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt)&
2911 &                                         -expipr(2)*radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt)
2912                                          radsumnfft1(2,klmn)=expipr(1)*radsum1(itypat)%value(2,klmn,i1,j2,i3,jkpt)&
2913 &                                         +expipr(2)*radsum1(itypat)%value(1,klmn,i1,j2,i3,jkpt)
2914                                          radsumnfft2(1,klmn)=expipr(1)*radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt)&
2915 &                                         -expipr(2)*radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt)
2916                                          radsumnfft2(2,klmn)=expipr(1)*radsum2(itypat)%value(2,klmn,i1,j2,i3,jkpt)&
2917 &                                         +expipr(2)*radsum2(itypat)%value(1,klmn,i1,j2,i3,jkpt)
2918                                          radsumnfft3(1,klmn)=expipr(1)*radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt)&
2919 &                                         -expipr(2)*radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt)
2920                                          radsumnfft3(2,klmn)=expipr(1)*radsum3(itypat)%value(2,klmn,i1,j2,i3,jkpt)&
2921 &                                         +expipr(2)*radsum3(itypat)%value(1,klmn,i1,j2,i3,jkpt)
2922                                        end do
2923 
2924 !                                      Loop on ij states
2925                                        do ilmn = 1, lmn_size
2926                                          i0lmn = ilmn*(ilmn-1)/2
2927                                          do jlmn = 1, lmn_size
2928                                            klmn = i0lmn+jlmn
2929                                            if (jlmn>ilmn) then
2930                                              i0lmn=jlmn*(jlmn-1)/2; klmn=i0lmn+ilmn
2931                                            end if
2932 !                                          Transform 3-dimentional radsum to 1-dimentional radsumnfft
2933                                            cp(:)=cprj_k(iatom,ib_cprj)%cp(:,ilmn)
2934                                            cp_pos(:)=cprj_pos(iatom,1)%cp(:,jlmn)
2935                                            cp11= cp(1)*cp_pos(1)
2936                                            cp22= cp(2)*cp_pos(2)*wf_fact
2937                                            cp12= cp(1)*cp_pos(2)*wf_fact
2938                                            cp21= cp(2)*cp_pos(1)
2939                                            cpr=cp11-cp22 ; cpi=cp12+cp21
2940                                            rho_contrib_paw1(1,indx) = rho_contrib_paw1(1,indx) &
2941 &                                           + four_pi*(cpr*radsumnfft1(1,klmn)-cpi*radsumnfft1(2,klmn))
2942                                            rho_contrib_paw1(2,indx) = rho_contrib_paw1(2,indx) &
2943 &                                           + four_pi*(cpr*radsumnfft1(2,klmn)+cpi*radsumnfft1(1,klmn))
2944                                            rho_contrib_paw2(1,indx) = rho_contrib_paw2(1,indx) &
2945 &                                           + four_pi*(cpr*radsumnfft2(1,klmn)-cpi*radsumnfft2(2,klmn))
2946                                            rho_contrib_paw2(2,indx) = rho_contrib_paw2(2,indx) &
2947 &                                           + four_pi*(cpr*radsumnfft2(2,klmn)+cpi*radsumnfft2(1,klmn))
2948                                            rho_contrib_paw3(1,indx) = rho_contrib_paw3(1,indx) &
2949 &                                           + four_pi*(cpr*radsumnfft3(1,klmn)-cpi*radsumnfft3(2,klmn))
2950                                            rho_contrib_paw3(2,indx) = rho_contrib_paw3(2,indx) &
2951 &                                           + four_pi*(cpr*radsumnfft3(2,klmn)+cpi*radsumnfft3(1,klmn))
2952                                          end do ! end loop over jlmn
2953                                        end do ! end loop over ilmn
2954 
2955                                      end do !end loop over atoms
2956 
2957                                      ABI_DEALLOCATE(radsumnfft1)
2958                                      ABI_DEALLOCATE(radsumnfft2)
2959                                      ABI_DEALLOCATE(radsumnfft3)
2960                                    end do !end loop over atom types
2961 
2962                                    rho_moment_v1(indx,jkpt) = rho_moment_v1(indx,jkpt) &
2963 &                                   +occ_el*occ_pos &
2964 &                                   *(gammastate*(rho_contrib_g(1,indx)**2+rho_contrib_g(2,indx)**2) &
2965 &                                   +rate_paw*(rho_contrib_paw1(1,indx)**2+rho_contrib_paw1(2,indx)**2 &
2966 &                                   -rho_contrib_paw2(1,indx)**2-rho_contrib_paw2(2,indx)**2))
2967                                    rho_moment_v2(indx,jkpt) = rho_moment_v2(indx,jkpt) &
2968 &                                   +occ_el*occ_pos*gammastate &
2969 &                                   *((rho_contrib_g(1,indx)+rho_contrib_paw3(1,indx))**2+&
2970 &                                   (rho_contrib_g(2,indx)+rho_contrib_paw3(2,indx))**2)
2971 
2972                                  end do ! end loop over i1
2973 
2974                                end if ! end loop over i2
2975                              end do
2976                            end do ! end loop over i3
2977 
2978                          end if ! PAW
2979 
2980 !                      ================================================================
2981 !                      End loops on electronic bands
2982 
2983                        end if ! occ>1.e-8
2984                      end do ! ibpp
2985                    end if ! occ_block>1.e-8
2986                  end do ! iblock
2987 
2988 !                End loops over k points and spins (electrons)
2989                  icg = icg + npw_k*my_nspinor*nband_k
2990                  ibg = ibg + my_nspinor*nband_cprj_k
2991                  ikg = ikg + npw_k
2992 
2993                  ABI_DEALLOCATE(cwaver)
2994                  ABI_DEALLOCATE(occ_k)
2995                  if (mpi_enreg%paral_kgb==0) then
2996                    ABI_DEALLOCATE(kg_k)
2997                    ABI_DEALLOCATE(gbound)
2998                  else
2999                    nullify(bandfft_kpt_el)
3000                  end if
3001                  if (dtset%usepaw==1) then
3002                    call pawcprj_free(cprj_k)
3003                    ABI_DATATYPE_DEALLOCATE(cprj_k)
3004                  end if
3005 
3006                end if ! mykpt
3007                bdtot_index=bdtot_index+nband_k
3008              end do ! ikpt
3009            end do ! isppol
3010 
3011 !          ================================================================
3012 !          End loops on positronic bands
3013 
3014          end if ! occ>1.e-8
3015        end do ! ibpp_pos
3016      end if ! occ(block)>1.e-8
3017    end do ! iblock_pos
3018 
3019 !  End loop over k points (positron)
3020    if (mykpt_pos) then
3021      icg_pos = icg_pos + npw_k_pos*my_nspinor*nband_k_pos
3022      ibg_pos = ibg_pos + my_nspinor*nband_cprj_k_pos
3023      ikg_pos = ikg_pos + npw_k_pos
3024    end if
3025    bdtot_index_pos=bdtot_index_pos+nband_k_pos
3026 
3027    ABI_DEALLOCATE(cwaver_pos)
3028    ABI_DEALLOCATE(cwaver_pos_block)
3029    ABI_DEALLOCATE(cg_k_pos)
3030    ABI_DEALLOCATE(occ_k_pos)
3031    if (mpi_enreg%paral_kgb==0) then
3032      ABI_DEALLOCATE(kg_k_pos)
3033      ABI_DEALLOCATE(gbound_pos)
3034    else if (mykpt_pos) then
3035      nullify(bandfft_kpt_pos)
3036    else
3037      call bandfft_kpt_destroy(bandfft_kpt_pos)
3038      ABI_DATATYPE_DEALLOCATE(bandfft_kpt_pos)
3039    end if
3040    if (dtset%usepaw==1) then
3041      call pawcprj_free(cprj_pos)
3042      ABI_DATATYPE_DEALLOCATE(cprj_pos)
3043      call pawcprj_free(cprj_k_pos)
3044      ABI_DATATYPE_DEALLOCATE(cprj_k_pos)
3045    end if
3046 
3047  end do ! ikpt_pos
3048 
3049 !================================================================
3050 !Final computations and printing
3051 
3052 !In case of parallelism, sum over the communicator(s)
3053  if (nproc_band>1) then
3054    ABI_ALLOCATE(mpibuf,(3*nfft,my_ngrid))
3055    do jkpt=1,my_ngrid
3056      mpibuf(       1:  nfft,jkpt)=rho_moment_v1(1:nfft,jkpt)
3057      mpibuf(  nfft+1:2*nfft,jkpt)=rho_moment_v2(1:nfft,jkpt)
3058      mpibuf(2*nfft+1:3*nfft,jkpt)=rho_pw       (1:nfft,jkpt)
3059    end do
3060    call xmpi_sum(mpibuf,mpi_enreg%comm_band,ierr)
3061    do jkpt=1,my_ngrid
3062      rho_moment_v1(1:nfft,jkpt)=mpibuf(       1:  nfft,jkpt)
3063      rho_moment_v2(1:nfft,jkpt)=mpibuf(  nfft+1:2*nfft,jkpt)
3064      rho_pw(1:nfft,jkpt)       =mpibuf(2*nfft+1:3*nfft,jkpt)
3065    end do
3066    ABI_DEALLOCATE(mpibuf)
3067  end if
3068  if (dtset%usepaw==1) then
3069    call xmpi_sum(rho_moment_core,mpi_enreg%comm_band,ierr)
3070  end if
3071 
3072 !Add valence and core contributions
3073  if (dtset%usepaw==1) then
3074    if (dtset%nsppol==2.and.my_nsppol==1) rho_moment_core(:,:)=half*rho_moment_core(:,:)
3075    rho_moment_v1(:,:)=rho_moment_v1(:,:)+rho_moment_core(:,:)
3076    rho_moment_v2(:,:)=rho_moment_v2(:,:)+rho_moment_core(:,:)
3077  end if
3078 
3079  units_=pi*(one/InvFineStruct)**3/Time_Sec/1.e12_dp/electronpositron%posocc
3080  scale_=(two_pi**2)/(Crystal%ucvol**two_thirds)
3081 
3082 !Integrate rho_moment over p
3083  buf(1)=sum(rho_moment_v1(1:nfft,1:my_ngrid))
3084  buf(2)=sum(rho_moment_v2(1:nfft,1:my_ngrid))
3085  buf(3)=sum(rho_moment_core(1:nfft,1:my_ngrid))
3086  buf(4)=sum(rho_pw(1:nfft,1:my_ngrid))
3087  call xmpi_sum(buf,mpi_enreg%comm_kpt,ierr)
3088  call xmpi_sum(buf,mpi_enreg%comm_fft,ierr)
3089  lambda_v1=buf(1)*units_/Crystal%ucvol/nkpt
3090  lambda_v2=buf(2)*units_/Crystal%ucvol/nkpt
3091  lambda_core=buf(3)*units_/Crystal%ucvol/nkpt
3092  lambda_pw=buf(4)*units_/Crystal%ucvol/nkpt
3093 
3094 !Write result in _DOPPLER file
3095 !Requires MPI-IO if nproc_fft>1
3096  if (me_band==0) then
3097    if (me_kpt==0) then
3098      filename_dop=trim(dtfil%filnam_ds(4))//'_DOPPLER'
3099      vec=sqrt(dot_product(Crystal%gprimd(:,3),Crystal%gprimd(:,3)))
3100      ABI_ALLOCATE(pcart_k,(3,nfft))
3101      ABI_ALLOCATE(rho_moment_k,(nfft))
3102      if (dtset%nsppol==2) then
3103        ABI_ALLOCATE(rho_moment_k2,(nfft))
3104      end if
3105      if (accessfil==IO_MODE_FORTRAN) then  ! >>>>> Fortran access
3106 !      Open file and write first line
3107        ierr=open_file(filename_dop,msg,newunit=unit_doppler,form='unformatted')
3108        write(unit_doppler) nfft,nkpt,Crystal%ucvol,Crystal%rprimd(:,:)
3109      else                                 ! >>>>> MPI-IO access
3110        unit_doppler=get_unit()
3111 !      Open file and write first line
3112        call WffOpen(IO_MODE_MPI,mpi_enreg%comm_fft,filename_dop,ierr,wff,0,me_fft,unit_doppler)
3113        if (me_fft==0) then
3114          call xderiveWRecInit(wff,ierr)
3115          call xderiveWrite(wff,n1*n2*n3,ierr)
3116          call xderiveWrite(wff,nkpt,ierr)
3117          call xderiveWrite(wff,Crystal%ucvol,ierr)
3118          call xderiveWrite(wff,Crystal%rprimd(:,:),ierr)
3119          call xderiveWRecEnd(wff,ierr)
3120        else
3121          call xmoveOff(wff,n_int=2,n_dp=10,n_mark=2)
3122        end if
3123 !      Store table of FFT points treated by current proc
3124        ABI_ALLOCATE(my_ffttab,(nfft))
3125        my_ffttab=0
3126        do i3=1,n3
3127          do i2=1,n2
3128            if (me_fft==fftn2_distrib(i2)) then
3129              indx0=n1*(n2*(i3-1)+(i2-1))
3130              indx=n1*(my_n2*(i3-1)+(ffti2_local(i2)-1))
3131              my_ffttab(indx+1:indx+n1)=(/(indx0+ii,ii=1,n1)/)
3132            end if
3133          end do
3134        end do
3135        ABI_ALLOCATE(mpibuf,(1,nfft))
3136      end if
3137    end if
3138 
3139    jkpt=0
3140    do ikpt=1,nkpt
3141      if (nproc_kpt==1) then
3142        rho_moment_k(1:nfft)=rho_moment_v2(1:nfft,ikpt)
3143      else
3144        if (my_gridtab(ikpt)/=0) jkpt=jkpt+1
3145        if (me_kpt==0) then
3146          if (my_gridtab(ikpt)==0) then
3147            tag=ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,1)
3148            call xmpi_recv(rho_moment_k,iproc,tag,mpi_enreg%comm_kpt,ierr)
3149            if (dtset%nsppol==2) then
3150              tag=2*ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,2)
3151              call xmpi_recv(rho_moment_k2,iproc,tag,mpi_enreg%comm_kpt,ierr)
3152              rho_moment_k(1:nfft)=rho_moment_k(1:nfft)+rho_moment_k2(1:nfft)
3153            end if
3154          else if (any(mpi_enreg%my_isppoltab(:)==1)) then
3155            rho_moment_k(1:nfft)=rho_moment_v2(1:nfft,jkpt)
3156            if (dtset%nsppol==2) then
3157              ii=2;if (mpi_enreg%my_isppoltab(2)==1) ii=1
3158              tag=ii*ikpt;iproc=mpi_enreg%proc_distrb(ikpt,1,ii)
3159              call xmpi_recv(rho_moment_k2,iproc,tag,mpi_enreg%comm_kpt,ierr)
3160              rho_moment_k(1:nfft)=rho_moment_k(1:nfft)+rho_moment_k2(1:nfft)
3161            end if
3162          end if
3163        else if (my_gridtab(ikpt)/=0) then
3164          if (mpi_enreg%my_isppoltab(1)==1) then
3165            tag=ikpt
3166            call xmpi_send(rho_moment_v2(1:nfft,jkpt),0,tag,mpi_enreg%comm_kpt,ierr)
3167          end if
3168          if (dtset%nsppol==2.and.mpi_enreg%my_isppoltab(2)==1) then
3169            tag=2*ikpt
3170            call xmpi_send(rho_moment_v2(1:nfft,jkpt),0,tag,mpi_enreg%comm_kpt,ierr)
3171          end if
3172        end if
3173      end if ! nproc_kpt>1
3174      if (me_kpt==0) then
3175        indx=0
3176        do i3=1,n3
3177          ig3=i3-(i3/id3)*n3-1
3178          do i2=1,n2
3179            if (me_fft/=fftn2_distrib(i2)) cycle
3180            ig2=i2-(i2/id2)*n2-1
3181            do i1=1,n1
3182              ig1=i1-(i1/id1)*n1-1
3183              indx=indx+1
3184              pcart_k(:,indx)=Crystal%gprimd(:,1)*real(ig1+dtset%kpt(1,ikpt))+&
3185 &             Crystal%gprimd(:,2)*real(ig2+dtset%kpt(2,ikpt))+&
3186 &             Crystal%gprimd(:,3)*real(ig3+dtset%kpt(3,ikpt))
3187            end do
3188          end do
3189        end do
3190        if (accessfil==IO_MODE_FORTRAN) then
3191          write(unit_doppler) pcart_k(1:3,1:nfft),rho_moment_k(1:nfft)
3192        else
3193          mpibuf(1,1:nfft)=rho_moment_k(1:nfft)
3194          call xderiveWRecInit(wff,ierr)
3195          call xderiveWrite(wff,pcart_k,3,nfft,mpi_enreg%comm_fft,my_ffttab,ierr)
3196          call xderiveWrite(wff,mpibuf ,1,nfft,mpi_enreg%comm_fft,my_ffttab,ierr)
3197          call xderiveWRecEnd(wff,ierr)
3198        end if
3199      end if
3200    end do
3201    if (me_kpt==0) then
3202      ABI_DEALLOCATE(pcart_k)
3203      ABI_DEALLOCATE(rho_moment_k)
3204      if (dtset%nsppol==2) then
3205        ABI_DEALLOCATE(rho_moment_k2)
3206      end if
3207      if (accessfil==IO_MODE_FORTRAN) then
3208        ierr=close_unit(unit_doppler,msg)
3209      else
3210        call WffClose(wff,ierr)
3211        ABI_DEALLOCATE(my_ffttab)
3212        ABI_DEALLOCATE(mpibuf)
3213      end if
3214    end if
3215  end if ! me_band==0
3216 
3217 !Write results
3218  write(msg,'(7a)') &
3219 & ' Computation of electron-positron pairs momentum distribution completed.',ch10,&
3220 & '-File ',trim(filename_dop),' has been created.',ch10,&
3221 & '-Use ~abinit/scripts/post_processing/posdopspectra.F90 to process it.'
3222  call wrtout(ab_out,msg,'COLL')
3223  call wrtout(std_out,msg,'COLL')
3224  msg=' Some e-p annihilation rates (ns-1) obtained by integration of e-p pairs momentum distribution:'
3225  call wrtout(std_out,msg,'COLL')
3226  write(msg,'(a,es22.12,3(2a,es22.12))') &
3227 & '   Lambda (from module of sum of PAW contrib.)  = ',lambda_v2*1000._dp,ch10,&
3228 & '     = lambda_core: ',lambda_core*1000._dp,ch10,&
3229 & '      +lambda_pw  : ',lambda_pw*1000._dp,ch10,&
3230 & '      +lambda_paw : ',(lambda_v2-lambda_core-lambda_pw)*1000._dp
3231  call wrtout(std_out,msg,'COLL')
3232  write(msg,'(4(a,es22.12,a))') &
3233 & '   Lambda (from sum of modules of PAW contrib.) = ',lambda_v1*1000._dp,ch10,&
3234 & '     = lambda_core: ',lambda_core*1000._dp,ch10,&
3235 & '      +lambda_pw  : ',lambda_pw*1000._dp,ch10,&
3236 & '      +lambda_paw : ',(lambda_v1-lambda_core-lambda_pw)*1000._dp,ch10
3237  call wrtout(std_out,msg,'COLL')
3238  write(msg,'(4a,es22.12,2a)') ch10,&
3239 & ' Annihilation rate obtained from integration of e-p pairs momentum distribution:',ch10,&
3240 & '   lambda=',lambda_v2*1000._dp,' ns-1',ch10
3241  call wrtout(ab_out,msg,'COLL')
3242 
3243 !Deallocate remaining memory
3244  ABI_DEALLOCATE(my_gridtab)
3245  ABI_DEALLOCATE(rho_pw)
3246  ABI_DEALLOCATE(rho_moment_v1)
3247  ABI_DEALLOCATE(rho_moment_v2)
3248  ABI_DEALLOCATE(rho_moment_core)
3249  ABI_DEALLOCATE(rho_contrib)
3250  ABI_DEALLOCATE(rho_contrib_g)
3251  ABI_DEALLOCATE(rho_contrib_paw1)
3252  ABI_DEALLOCATE(rho_contrib_paw2)
3253  ABI_DEALLOCATE(rho_contrib_paw3)
3254  if (state_dependent) then
3255    call unset_mpi_enreg_fft(mpi_enreg_seq)
3256    call destroy_mpi_enreg(mpi_enreg_seq)
3257    ABI_DEALLOCATE(rhor_dop_el)
3258    if (dtset%usepaw==1) then
3259      call pawrhoij_free(pawrhoij_dop_el)
3260      ABI_DATATYPE_DEALLOCATE(pawrhoij_dop_el)
3261      if (mpi_enreg%my_natom<dtset%natom) then
3262        call pawrhoij_free(pawrhoij_all)
3263        call pawrhoij_free(pawrhoij_ep_all)
3264        ABI_DATATYPE_DEALLOCATE(pawrhoij_all)
3265        ABI_DATATYPE_DEALLOCATE(pawrhoij_ep_all)
3266      end if
3267    end if
3268  end if
3269 
3270  ABI_DEALLOCATE(gamma)
3271 
3272  if (dtset%usepaw==1.and.(.not.include_nhat_in_gamma)) then
3273    ABI_DEALLOCATE(rhor_)
3274    ABI_DEALLOCATE(rhor_ep_)
3275  end if
3276 
3277  if (dtset%usepaw==1) then
3278    ABI_DEALLOCATE(nphicor)
3279    ABI_DEALLOCATE(lmncmax)
3280    do itypat=1,dtset%ntypat
3281      if (allocated(phicor(itypat)%value)) then
3282        ABI_DEALLOCATE(phicor(itypat)%value)
3283      end if
3284      if (allocated(indlmncor(itypat)%value)) then
3285        ABI_DEALLOCATE(indlmncor(itypat)%value)
3286      end if
3287      if (allocated(radsumc(itypat)%value)) then
3288        ABI_DEALLOCATE(radsumc(itypat)%value)
3289      end if
3290      if (allocated(radsum1(itypat)%value)) then
3291        ABI_DEALLOCATE(radsum1(itypat)%value)
3292      end if
3293      if (allocated(radsum2(itypat)%value)) then
3294        ABI_DEALLOCATE(radsum2(itypat)%value)
3295      end if
3296      if (allocated(radsum3(itypat)%value)) then
3297        ABI_DEALLOCATE(radsum3(itypat)%value)
3298      end if
3299    end do
3300    do iatom=1,dtset%natom
3301      if (allocated(gammastate_c(iatom)%value)) then
3302        ABI_DEALLOCATE(gammastate_c(iatom)%value)
3303      end if
3304    end do
3305    ABI_DATATYPE_DEALLOCATE(phicor)
3306    ABI_DATATYPE_DEALLOCATE(indlmncor)
3307    ABI_DATATYPE_DEALLOCATE(radsumc)
3308    ABI_DATATYPE_DEALLOCATE(radsum1)
3309    ABI_DATATYPE_DEALLOCATE(radsum2)
3310    ABI_DATATYPE_DEALLOCATE(radsum3)
3311    ABI_DATATYPE_DEALLOCATE(gammastate_c)
3312  end if
3313 
3314  DBG_EXIT("COLL")
3315 
3316 end subroutine posdoppler

ABINIT/poslifetime [ Functions ]

[ Top ] [ 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

PARENTS

      outscfcv,posdoppler

CHILDREN

      gammapositron,gammapositron_fft,mkdenpos,nderiv_gen,pawdensities
      pawxcsum,simp_gen,wrtout,xmpi_sum

SOURCE

 915 subroutine poslifetime(dtset,electronpositron,gprimd,my_natom,mpi_enreg,n3xccc,nfft,ngfft,nhat,&
 916 &                      option,pawang,pawrad,pawrhoij,pawtab,rate,rate_paw,rhor,ucvol,xccc3d,&
 917 &                      rhor_dop_el,pawrhoij_dop_el,pawrhoij_ep) ! optional arguments
 918 
 919 
 920 !This section has been created automatically by the script Abilint (TD).
 921 !Do not modify the following lines by hand.
 922 #undef ABI_FUNC
 923 #define ABI_FUNC 'poslifetime'
 924 !End of the abilint section
 925 
 926  implicit none
 927 
 928 !Arguments ------------------------------------
 929 !scalars
 930  integer,intent(in) :: my_natom,n3xccc,nfft,option
 931  real(dp),intent(in) :: ucvol
 932  real(dp),intent(out) :: rate,rate_paw
 933  type(dataset_type), intent(in) :: dtset
 934  type(electronpositron_type),pointer :: electronpositron
 935  type(MPI_type),intent(in) :: mpi_enreg
 936  type(pawang_type), intent(in) :: pawang
 937 !arrays
 938  integer,intent(in) :: ngfft(18)
 939  real(dp),intent(in) :: gprimd(3,3),nhat(nfft,dtset%nspden*dtset%usepaw),xccc3d(n3xccc)
 940  real(dp),intent(in),target :: rhor(nfft,dtset%nspden)
 941  real(dp),optional,intent(in) :: rhor_dop_el(nfft)
 942  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw)
 943  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*dtset%usepaw)
 944  type(pawrhoij_type),optional,intent(in) :: pawrhoij_dop_el(my_natom*dtset%usepaw)
 945  type(pawrhoij_type),optional,target,intent(in) :: pawrhoij_ep(my_natom*dtset%usepaw)
 946  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
 947 
 948 !Local variables-------------------------------
 949 !scalars
 950  integer :: cplex,iatom,ierr,ifft,igam,ii,ilm,ilm1,ilm2,iloop,ipt,ir,isel
 951  integer :: itypat,iwarn,iwarnj,iwarnp,lm_size,lmn2_size,mesh_size
 952  integer :: nfftot,ngamma,ngr,ngrad,nspden_ep,opt_dens,usecore
 953  logical,parameter :: include_nhat_in_gamma=.false.
 954  real(dp),parameter :: delta=1.d-4
 955  real(dp) :: fact,fact2,intg
 956  real(dp) :: lambda_core    ,lambda_core_ipm    ,lambda    ,lambda_ipm
 957  real(dp) :: lambda_core_paw,lambda_core_paw_ipm,lambda_paw,lambda_paw_ipm
 958  real(dp) :: lifetime,lifetime_ipm,nbec,nbev,nbp,rdum,sqfpi,units
 959  character(len=500) :: msg
 960 !arrays
 961  integer,allocatable :: igamma(:)
 962  logical,allocatable :: lmselect(:),lmselect_ep(:),lmselect_dum(:)
 963  real(dp) :: mpibuf(4)
 964  real(dp),parameter :: qphon(3)=(/zero,zero,zero/),lsign(2)=(/one,-one/)
 965  real(dp),allocatable :: d1gam(:,:,:),d2gam(:,:,:),ff(:),gam_(:,:,:),gamma(:,:),gammam(:,:,:),gg(:,:)
 966  real(dp),allocatable :: grhocore2(:),grhocor2_(:),grhoe2(:),grho2_(:)
 967  real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:),nhat1_j(:,:,:)
 968  real(dp),allocatable :: rho_(:),rho_ep_(:),rho1(:,:,:),rho1_ep(:,:,:),rho1_j(:,:,:)
 969  real(dp),allocatable :: rhoarr1(:),rhoarr1_ep(:),rhoarr1_j(:),rhoarr2(:)
 970  real(dp),allocatable :: rhocore(:),rhocor_(:),rhoe(:,:),rhop(:,:),rhor_dop_el_(:)
 971  real(dp),allocatable :: rhosph(:),rhosph_ep(:),rhosph_j(:),rhotot(:,:),rhotot_ep(:,:)
 972  real(dp),allocatable :: rhotot_j(:,:),trho1(:,:,:),trho1_ep(:,:,:),trho1_j(:,:,:)
 973  real(dp),allocatable :: v1sum(:,:),v2sum(:,:,:)
 974  real(dp),pointer :: rhor_(:,:),rhor_ep_(:,:)
 975  type(pawrhoij_type),pointer :: pawrhoij_ep_(:)
 976 
 977 ! *************************************************************************
 978 
 979  DBG_ENTER("COLL")
 980 
 981 !Tests for developers
 982  if (.not.associated(electronpositron)) then
 983    msg='electronpositron variable must be associated!'
 984    MSG_BUG(msg)
 985  end if
 986  if (option/=1) then
 987    if ((.not.present(rhor_dop_el)).or.(.not.present(pawrhoij_dop_el))) then
 988      msg='when option/=1, rhor_dop_el and pawrhoij_dop_el must be present!'
 989      MSG_BUG(msg)
 990    end if
 991  end if
 992 
 993 !Constants
 994  fact=0.0
 995  cplex=1;nspden_ep=1
 996  usecore=n3xccc/nfft
 997  nfftot=ngfft(1)*ngfft(2)*ngfft(3)
 998  ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
 999  iwarn=0;iwarnj=0;iwarnp=1
1000  sqfpi=sqrt(four_pi)
1001 
1002 !Compatibility tests
1003  if (electronpositron%particle==EP_NOTHING) then
1004    msg='Not valid for electronpositron%particle=NOTHING!'
1005    MSG_BUG(msg)
1006  end if
1007  if (electronpositron%nfft/=nfft) then
1008    msg='nfft/=electronpositron%nfft!'
1009    MSG_BUG(msg)
1010  end if
1011  if (dtset%usepaw==1) then
1012    if(dtset%pawxcdev==0.and.ngrad==2) then
1013      msg='GGA is not implemented for pawxcdev=0 (use dtset%pawxcdev/=0)!'
1014      MSG_BUG(msg)
1015    end if
1016  end if
1017 
1018 !Select type(s) of enhancement factor
1019  if ((electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3).and.option==1) then
1020    ngamma=2
1021    ABI_ALLOCATE(igamma,(ngamma))
1022    igamma(1)=1;igamma(2)=2
1023  else
1024    ngamma=1
1025    ABI_ALLOCATE(igamma,(ngamma))
1026    if (electronpositron%ixcpositron==-1) igamma(1)=0
1027    if (electronpositron%ixcpositron== 2) igamma(1)=4
1028    if (electronpositron%ixcpositron==11.or.electronpositron%ixcpositron==31) igamma(1)=3
1029    if (electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3) igamma(1)=2
1030  end if
1031 
1032 !Select density according to nhat choice
1033  if (dtset%usepaw==0.or.include_nhat_in_gamma) then
1034    rhor_ => rhor
1035    rhor_ep_ => electronpositron%rhor_ep
1036  else
1037    ABI_ALLOCATE(rhor_,(nfft,dtset%nspden))
1038    ABI_ALLOCATE(rhor_ep_,(nfft,dtset%nspden))
1039    rhor_=rhor-nhat
1040    rhor_ep_=electronpositron%rhor_ep-electronpositron%nhat_ep
1041  end if
1042 
1043 !Eventually overwrite electronpositron%pawrhoij_ep
1044  if (present(pawrhoij_ep)) then
1045    pawrhoij_ep_ => pawrhoij_ep
1046  else
1047    pawrhoij_ep_ => electronpositron%pawrhoij_ep
1048  end if
1049 
1050 !Loop on different enhancement factors
1051  do igam=1,ngamma
1052 
1053 !  Compute electron-positron annihilation rate using pseudo densities (plane waves)
1054 !  ----------------------------------------------------------------------------------------
1055 
1056 !  Select the densities and make them positive
1057    ABI_ALLOCATE(rhoe,(nfft,nspden_ep))
1058    ABI_ALLOCATE(rhop,(nfft,nspden_ep))
1059    if (electronpositron%particle==EP_ELECTRON) then
1060      rhoe(:,1)=rhor_ep_(:,1);rhop(:,1)=rhor_(:,1)
1061    else if (electronpositron%particle==EP_POSITRON) then
1062      rhoe(:,1)=rhor_(:,1);rhop(:,1)=rhor_ep_(:,1)
1063    end if
1064    call mkdenpos(iwarn ,nfft,nspden_ep,1,rhoe,dtset%xc_denpos)
1065    call mkdenpos(iwarnp,nfft,nspden_ep,1,rhop,dtset%xc_denpos)
1066    if (option/=1) then
1067      ABI_ALLOCATE(rhor_dop_el_,(nfft))
1068      rhor_dop_el_(:)=rhor_dop_el(:)
1069      call mkdenpos(iwarnp,nfft,1,1,rhor_dop_el_,dtset%xc_denpos)
1070    end if
1071 
1072 !  Compute enhancement factor at each FFT grid point
1073 !  gamma(:,1): using total   electronic density
1074 !  gamma(:,2): using valence electronic density
1075    ABI_ALLOCATE(gamma,(nfft,2))
1076    if (option==1.or.option==2) then
1077      call gammapositron_fft(electronpositron,gamma,gprimd,igamma(igam),mpi_enreg,&
1078 &     n3xccc,nfft,ngfft,rhoe,rhop,xccc3d)
1079    else
1080      gamma=one
1081    end if
1082 
1083 !  Compute positron annihilation rates
1084    lambda     =zero;lambda_ipm     =zero
1085    lambda_core=zero;lambda_core_ipm=zero
1086    if (option==1) then
1087      do ifft=1,nfft
1088        lambda    =lambda    +rhop(ifft,1)*rhoe(ifft,1)*gamma(ifft,1)
1089        lambda_ipm=lambda_ipm+rhop(ifft,1)*rhoe(ifft,1)*gamma(ifft,2)
1090      end do
1091    else
1092      do ifft=1,nfft
1093        lambda    =lambda    +rhop(ifft,1)*rhor_dop_el_(ifft)*gamma(ifft,1)
1094        lambda_ipm=lambda_ipm+rhop(ifft,1)*rhor_dop_el_(ifft)*gamma(ifft,2)
1095      end do
1096    end if
1097    if (usecore==1) then
1098      do ifft=1,nfft
1099        lambda_core    =lambda_core    +rhop(ifft,1)*xccc3d(ifft)*gamma(ifft,1)
1100        lambda_core_ipm=lambda_core_ipm+rhop(ifft,1)*xccc3d(ifft)
1101      end do
1102    end if
1103    lambda         =lambda         *ucvol/dble(nfftot)
1104    lambda_ipm     =lambda_ipm     *ucvol/dble(nfftot)
1105    lambda_core    =lambda_core    *ucvol/dble(nfftot)
1106    lambda_core_ipm=lambda_core_ipm*ucvol/dble(nfftot)
1107    ABI_DEALLOCATE(gamma)
1108    ABI_DEALLOCATE(rhoe)
1109    ABI_DEALLOCATE(rhop)
1110    if (option/=1) then
1111      ABI_DEALLOCATE(rhor_dop_el_)
1112    end if
1113 !  NC pseudopotential: check electrons/positron number
1114    if (dtset%usepaw==0.and.igam==ngamma) then
1115      nbec=zero;nbev=zero;nbp=zero
1116      if (electronpositron%particle==EP_ELECTRON) then
1117        do ifft=1,nfft
1118          nbec=nbec+xccc3d(ifft)
1119          nbev=nbev+electronpositron%rhor_ep(ifft,1)
1120          nbp =nbp +rhor(ifft,1)
1121        end do
1122      else
1123        do ifft=1,nfft
1124          nbec=nbec+xccc3d(ifft)
1125          nbev=nbev+rhor(ifft,1)
1126          nbp =nbp +electronpositron%rhor_ep(ifft,1)
1127        end do
1128      end if
1129      nbec=nbec*ucvol/dble(nfftot)
1130      nbev=nbev*ucvol/dble(nfftot)
1131      nbp =nbp *ucvol/dble(nfftot)
1132    end if
1133 
1134 !  MPI parallelization
1135    if(mpi_enreg%nproc_fft>1)then
1136      call xmpi_sum(lambda    ,mpi_enreg%comm_fft,ierr)
1137      call xmpi_sum(lambda_ipm,mpi_enreg%comm_fft,ierr)
1138      call xmpi_sum(lambda_core    ,mpi_enreg%comm_fft,ierr)
1139      call xmpi_sum(lambda_core_ipm,mpi_enreg%comm_fft,ierr)
1140      if (dtset%usepaw==0.and.igam==ngamma) then
1141        call xmpi_sum(nbec,mpi_enreg%comm_fft,ierr)
1142        call xmpi_sum(nbev,mpi_enreg%comm_fft,ierr)
1143        call xmpi_sum(nbp ,mpi_enreg%comm_fft,ierr)
1144      end if
1145    end if
1146 
1147 
1148 !  PAW: add on-site contributions to electron-positron annihilation rate
1149 !  ----------------------------------------------------------------------------------------
1150    if (dtset%usepaw==1) then
1151 
1152      lambda_paw     =zero;lambda_paw_ipm     =zero
1153      lambda_core_paw=zero;lambda_core_paw_ipm=zero
1154 
1155 !    Loop on atoms
1156      do iatom=1,my_natom
1157 
1158        itypat=pawrhoij(iatom)%itypat
1159        lmn2_size=pawtab(itypat)%lmn2_size
1160        mesh_size=pawtab(itypat)%mesh_size
1161        lm_size=pawtab(itypat)%lcut_size**2
1162        cplex=1
1163        ngr=0;if (ngrad==2) ngr=mesh_size
1164 
1165 !      Allocations of "on-site" densities
1166        ABI_ALLOCATE(rho1 ,(cplex*mesh_size,lm_size,nspden_ep))
1167        ABI_ALLOCATE(trho1,(cplex*mesh_size,lm_size,nspden_ep))
1168        ABI_ALLOCATE(rho1_ep ,(cplex*mesh_size,lm_size,nspden_ep))
1169        ABI_ALLOCATE(trho1_ep,(cplex*mesh_size,lm_size,nspden_ep))
1170        if (option/=1) then
1171          ABI_ALLOCATE(rho1_j ,(cplex*mesh_size,lm_size,nspden_ep))
1172          ABI_ALLOCATE(trho1_j,(cplex*mesh_size,lm_size,nspden_ep))
1173        else
1174          ABI_ALLOCATE(rho1_j ,(0,0,0))
1175          ABI_ALLOCATE(trho1_j,(0,0,0))
1176        end if
1177        if (include_nhat_in_gamma) then
1178          ABI_ALLOCATE(nhat1,(cplex*mesh_size,lm_size,nspden_ep))
1179          ABI_ALLOCATE(nhat1_ep,(cplex*mesh_size,lm_size,nspden_ep))
1180        else
1181          ABI_ALLOCATE(nhat1,(0,0,0))
1182          ABI_ALLOCATE(nhat1_ep,(0,0,0))
1183        end if
1184        if (include_nhat_in_gamma.and.option/=1) then
1185          ABI_ALLOCATE(nhat1_j,(cplex*mesh_size,lm_size,nspden_ep))
1186        else
1187          ABI_ALLOCATE(nhat1_j,(0,0,0))
1188        end if
1189        ABI_ALLOCATE(lmselect,(lm_size))
1190        ABI_ALLOCATE(lmselect_ep,(lm_size))
1191        ABI_ALLOCATE(lmselect_dum,(lm_size))
1192 
1193 !      Compute "on-site" densities (n1, ntild1, nhat1) for electron and positron =====
1194        lmselect(:)=.true.
1195        opt_dens=1;if (include_nhat_in_gamma) opt_dens=0
1196        call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1,nspden_ep,1,&
1197 &       0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij(iatom),&
1198 &       pawtab(itypat),rho1,trho1)
1199        lmselect_ep(:)=.true.
1200        call pawdensities(rdum,cplex,iatom,lmselect_ep,lmselect_dum,lm_size,nhat1_ep,nspden_ep,1,&
1201 &       0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_ep_(iatom),&
1202 &       pawtab(itypat),rho1_ep,trho1_ep)
1203 
1204 !      For state dependent scheme in Doppler                                       =====
1205 !      Compute "on-site" densities (n1, ntild1, nhat1) for a given electron state j=====
1206        if (option/=1) then
1207          opt_dens=1;if (include_nhat_in_gamma) opt_dens=0
1208          call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1_j,nspden_ep,1,&
1209 &         0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_dop_el(iatom),&
1210 &         pawtab(itypat),rho1_j,trho1_j)
1211        end if
1212 
1213 !      Compute contribution to annihilation rate:
1214 !      Loop: first step: compute all-electron contribution (from n^1, n_c)
1215 !      2nd   step: compute pseudo contribution (from tild_n^1, hat_n^1, tild_n_c)
1216        do iloop=1,2
1217          if (iloop==1) usecore=1
1218          if (iloop==2) usecore=pawtab(itypat)%usetcore
1219          ABI_ALLOCATE(rhocore,(mesh_size))
1220 
1221 !        First formalism: use densities on r,theta,phi
1222          if (dtset%pawxcdev==0) then
1223 
1224            ABI_ALLOCATE(gamma,(mesh_size,2))
1225            ABI_ALLOCATE(rhoarr1,(mesh_size))
1226            ABI_ALLOCATE(rhoarr1_ep,(mesh_size))
1227            if (option/=1) then
1228              ABI_ALLOCATE(rhoarr1_j,(mesh_size))
1229            end if
1230 !          Loop on the angular part
1231            do ipt=1,pawang%angl_size
1232 !            Build densities
1233              rhoarr1=zero;rhoarr1_ep=zero;rhocore=zero
1234              if (option/=1) rhoarr1_j=zero
1235              if (iloop==1) then
1236                do ilm=1,lm_size
1237                  if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+rho1(:,ilm,1)*pawang%ylmr(ilm,ipt)
1238                end do
1239                if (option/=1) then
1240                  do ilm=1,lm_size
1241                    if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1_j(:)+rho1_j(:,ilm,1)*pawang%ylmr(ilm,ipt)
1242                  end do
1243                end if
1244                do ilm=1,lm_size
1245                  if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+rho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt)
1246                end do
1247                if (usecore==1) rhocore(:)=pawtab(itypat)%coredens(:)
1248              else
1249                if (include_nhat_in_gamma) then
1250                  do ilm=1,lm_size
1251                    if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+(trho1(:,ilm,1)+nhat1(:,ilm,1))*pawang%ylmr(ilm,ipt)
1252                  end do
1253                  if (option/=1) then
1254                    do ilm=1,lm_size
1255                      if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1_j(:)+(trho1_j(:,ilm,1)+nhat1_j(:,ilm,1))*pawang%ylmr(ilm,ipt)
1256                    end do
1257                  end if
1258                  do ilm=1,lm_size
1259                    if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+(trho1_ep(:,ilm,1)+nhat1_ep(:,ilm,1))*pawang%ylmr(ilm,ipt)
1260                  end do
1261                else
1262                  do ilm=1,lm_size
1263                    if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+trho1(:,ilm,1)*pawang%ylmr(ilm,ipt)
1264                  end do
1265                  if (option/=1) then
1266                    do ilm=1,lm_size
1267                      if (lmselect(ilm)) rhoarr1_j(:)=rhoarr1(:)+trho1_j(:,ilm,1)*pawang%ylmr(ilm,ipt)
1268                    end do
1269                  end if
1270                  do ilm=1,lm_size
1271                    if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+trho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt)
1272                  end do
1273                end if
1274                if (usecore==1) rhocore(:)=pawtab(itypat)%tcoredens(:,1)
1275              end if
1276 !            Make the densities positive
1277              if (electronpositron%particle==EP_ELECTRON) then
1278                call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1   ,dtset%xc_denpos)
1279                call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos)
1280              else if (electronpositron%particle==EP_POSITRON) then
1281                call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1   ,dtset%xc_denpos)
1282                call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos)
1283                if (option/=1) then
1284                  call mkdenpos(iwarnj,mesh_size,1,1,rhoarr1_j,dtset%xc_denpos)
1285                end if
1286              end if
1287 !            Compute Gamma
1288              ABI_ALLOCATE(grhoe2,(ngr))
1289              ABI_ALLOCATE(grhocore2,(ngr))
1290              if (option==1.or.option==2) then
1291                if (electronpositron%particle==EP_ELECTRON) then
1292                  call gammapositron(gamma,grhocore2,grhoe2,igamma(igam),ngr,mesh_size,&
1293 &                 rhocore,rhoarr1_ep,rhoarr1,usecore)
1294                else if (electronpositron%particle==EP_POSITRON) then
1295                  call gammapositron(gamma,grhocore2,grhoe2,igamma(igam),ngr,mesh_size,&
1296 &                 rhocore,rhoarr1,rhoarr1_ep,usecore)
1297                end if
1298              else
1299                gamma(:,:)=one
1300              end if
1301              ABI_DEALLOCATE(grhoe2)
1302              ABI_DEALLOCATE(grhocore2)
1303 !            Compute contribution to annihilation rates
1304              ABI_ALLOCATE(ff,(mesh_size))
1305              if (option/=1) rhoarr1(:)=rhoarr1_j(:)
1306              do ii=1,4
1307                if (ii==1) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhoarr1_ep(1:mesh_size) &
1308 &               *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
1309                if (ii==2) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhoarr1_ep(1:mesh_size) &
1310 &               *gamma(1:mesh_size,2)*pawrad(itypat)%rad(1:mesh_size)**2
1311                if (electronpositron%particle==EP_ELECTRON) then
1312                  if (ii==3) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhocore(1:mesh_size) &
1313 &                 *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
1314                  if (ii==4) ff(1:mesh_size)=rhoarr1(1:mesh_size)*rhocore(1:mesh_size) &
1315 &                 *pawrad(itypat)%rad(1:mesh_size)**2
1316                else
1317                  if (ii==3) ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocore(1:mesh_size) &
1318 &                 *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
1319                  if (ii==4) ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocore(1:mesh_size) &
1320 &                 *pawrad(itypat)%rad(1:mesh_size)**2
1321                end if
1322                call simp_gen(intg,ff,pawrad(itypat))
1323                intg=intg*pawang%angwgth(ipt)*four_pi
1324                if (ii==1) lambda_paw         =lambda_paw         +lsign(iloop)*intg
1325                if (ii==2) lambda_paw_ipm     =lambda_paw_ipm     +lsign(iloop)*intg
1326                if (ii==3) lambda_core_paw    =lambda_core_paw    +lsign(iloop)*intg
1327                if (ii==4) lambda_core_paw_ipm=lambda_core_paw_ipm+lsign(iloop)*intg
1328              end do
1329              ABI_DEALLOCATE(ff)
1330            end do ! ipt
1331            ABI_DEALLOCATE(gamma)
1332            ABI_DEALLOCATE(rhoarr1)
1333            ABI_DEALLOCATE(rhoarr1_ep)
1334            if (option/=1) then
1335              ABI_DEALLOCATE(rhoarr1_j)
1336            end if
1337 
1338 !          Second formalism: use (l,m) moments for densities
1339          else if (dtset%pawxcdev/=0) then
1340 
1341 !          Build densities
1342            ABI_ALLOCATE(gammam,(mesh_size,2,lm_size))
1343            ABI_ALLOCATE(rhotot,(mesh_size,lm_size))
1344            ABI_ALLOCATE(rhotot_ep,(mesh_size,lm_size))
1345            ABI_ALLOCATE(rhosph,(mesh_size))
1346            ABI_ALLOCATE(rhosph_ep,(mesh_size))
1347            if (option/=1) then
1348              ABI_ALLOCATE(rhotot_j,(mesh_size,lm_size))
1349              ABI_ALLOCATE(rhosph_j,(mesh_size))
1350            end if
1351            if (usecore==0) rhocore(:)=zero
1352            if (iloop==1) then
1353              rhotot   (:,:)=rho1   (:,:,1)
1354              rhotot_ep(:,:)=rho1_ep(:,:,1)
1355              if (option/=1) rhotot_j (:,:)=rho1_j (:,:,1)
1356              if (usecore==1) rhocore(:)=pawtab(itypat)%coredens(:)
1357            else
1358              if (include_nhat_in_gamma) then
1359                rhotot   (:,:)=trho1   (:,:,1)+nhat1   (:,:,1)
1360                rhotot_ep(:,:)=trho1_ep(:,:,1)+nhat1_ep(:,:,1)
1361                if (option/=1) rhotot_j (:,:)=trho1_j (:,:,1)+nhat1_j (:,:,1)
1362              else
1363                rhotot   (:,:)=trho1   (:,:,1)
1364                rhotot_ep(:,:)=trho1_ep(:,:,1)
1365                if (option/=1) rhotot_j (:,:)=trho1_j (:,:,1)
1366              end if
1367              if (usecore==1) rhocore(:)=pawtab(itypat)%tcoredens(:,1)
1368            end if
1369            rhosph   (:)=rhotot   (:,1)/sqfpi
1370            rhosph_ep(:)=rhotot_ep(:,1)/sqfpi
1371            if (option/=1) rhosph_j (:)=rhotot_j(:,1)/sqfpi
1372 !          Make spherical densities positive
1373            if (electronpositron%particle==EP_ELECTRON) then
1374              call mkdenpos(iwarnp,mesh_size,1,1,rhosph   ,dtset%xc_denpos)
1375              call mkdenpos(iwarn ,mesh_size,1,1,rhosph_ep,dtset%xc_denpos)
1376            else if (electronpositron%particle==EP_POSITRON) then
1377              call mkdenpos(iwarn ,mesh_size,1,1,rhosph   ,dtset%xc_denpos)
1378              call mkdenpos(iwarnp,mesh_size,1,1,rhosph_ep,dtset%xc_denpos)
1379              if (option/=1) then
1380                call mkdenpos(iwarnp,mesh_size,1,1,rhosph_j,dtset%xc_denpos)
1381              end if
1382            end if
1383 !          Need gradients of electronic densities for GGA
1384            ABI_ALLOCATE(grhoe2,(ngr))
1385            ABI_ALLOCATE(grhocore2,(ngr))
1386            if (ngr>0) then
1387              if (electronpositron%particle==EP_ELECTRON) then
1388                call nderiv_gen(grhoe2,rhosph_ep,pawrad(itypat))
1389              else if (electronpositron%particle==EP_POSITRON) then
1390                call nderiv_gen(grhoe2,rhosph,pawrad(itypat))
1391              end if
1392              grhoe2(:)=grhoe2(:)**2
1393              if (usecore==1) then
1394                call nderiv_gen(grhocore2,rhocore,pawrad(itypat))
1395                grhocore2(:)=grhocore2(:)**2
1396              end if
1397            end if
1398 !          Compute Gamma for (rho-,rho+),
1399 !          (rho- +drho-,rho+), (rho- -drho-,rho+),
1400 !          (rho-,rho+ +drho+), (rho-,rho+ -drho+),
1401 !          (rho- +drho-,rho+ +drho+), (rho- -drho-,rho+ -drho+)
1402 !          Do a seven steps loop
1403            ABI_ALLOCATE(gam_,(mesh_size,2,7))
1404            ABI_ALLOCATE(rho_,(mesh_size))
1405            ABI_ALLOCATE(rho_ep_,(mesh_size))
1406            ABI_ALLOCATE(rhocor_,(mesh_size))
1407            ABI_ALLOCATE(grho2_,(ngr))
1408            ABI_ALLOCATE(grhocor2_,(ngr))
1409            do ii=1,7
1410 !            Apply delta to get perturbed densities
1411              rho_(:)=rhosph(:);rho_ep_(:)=rhosph_ep(:);if (usecore==1) rhocor_(:)=rhocore(:)
1412              if (ngr>0) grho2_(:)=grhoe2(:)
1413              if (ngr>0) grhocor2_(:)=grhocore2(:)
1414              if (ii==2.or.ii==4.or.ii==6) fact=(one+delta)
1415              if (ii==3.or.ii==5.or.ii==7) fact=(one-delta)
1416              fact2=fact**2
1417              if (ii==2.or.ii==3.or.ii==6.or.ii==7) then
1418                rho_(:)=fact*rho_(:)
1419                if (electronpositron%particle==EP_POSITRON) then
1420                  if (ngr>0) grho2_(:)=fact2*grho2_(:)
1421                  if (usecore==1)rhocor_(:)=fact*rhocor_(:)
1422                  if (ngr>0.and.usecore==1) grhocor2_(:)=fact2*grhocor2_(:)
1423                end if
1424              end if
1425              if (ii==4.or.ii==5.or.ii==6.or.ii==7) then
1426                rho_ep_(:)=fact*rho_ep_(:)
1427                if (electronpositron%particle==EP_ELECTRON) then
1428                  if (ngr>0) grho2_(:)=fact2*grho2_(:)
1429                  if (usecore==1)rhocor_(:)=fact*rhocor_(:)
1430                  if (ngr>0.and.usecore==1) grhocor2_(:)=fact2*grhocor2_(:)
1431                end if
1432              end if
1433 !            Compute gamma for these perturbed densities
1434              if (option==1.or.option==2) then
1435                if (electronpositron%particle==EP_ELECTRON) then
1436                  call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma(igam),ngr,mesh_size,rhocor_,rho_ep_,rho_,usecore)
1437                else if (electronpositron%particle==EP_POSITRON) then
1438                  call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma(igam),ngr,mesh_size,rhocor_,rho_,rho_ep_,usecore)
1439                end if
1440              else
1441                gam_(:,:,:)=one
1442              end if
1443            end do ! end loop ii=1,7
1444 
1445            ABI_DEALLOCATE(rhocor_)
1446            ABI_DEALLOCATE(grho2_)
1447            ABI_DEALLOCATE(grhocor2_)
1448            ABI_DEALLOCATE(grhoe2)
1449            ABI_DEALLOCATE(grhocore2)
1450            rho_   (:)=rhosph   (:);if (electronpositron%particle==EP_POSITRON.and.usecore==1) rho_   (:)=rho_   (:)+rhocore(:)
1451            rho_ep_(:)=rhosph_ep(:);if (electronpositron%particle==EP_ELECTRON.and.usecore==1) rho_ep_(:)=rho_ep_(:)+rhocore(:)
1452 !          Compute numerical first and second derivatives of Gamma
1453 !          d1gam(1) = dgam/drho+ (particle=ELECTRON), dgam/drho- (particle=POSITRON)
1454 !          d1gam(2) = dgam/drho- (particle=ELECTRON), dgam/drho+ (particle=POSITRON)
1455            ABI_ALLOCATE(d1gam,(mesh_size,2,2))
1456            d1gam(:,:,:)=zero
1457            do ir=1,mesh_size
1458              if (rho_     (ir)>tol14) d1gam(ir,1,1)=(gam_(ir,1,2)-gam_(ir,1,3))*half/(delta*rho_     (ir))
1459              if (rhosph   (ir)>tol14) d1gam(ir,2,1)=(gam_(ir,2,2)-gam_(ir,2,3))*half/(delta*rhosph   (ir))
1460              if (rho_ep_  (ir)>tol14) d1gam(ir,1,2)=(gam_(ir,1,4)-gam_(ir,1,5))*half/(delta*rho_ep_  (ir))
1461              if (rhosph_ep(ir)>tol14) d1gam(ir,2,2)=(gam_(ir,2,4)-gam_(ir,2,5))*half/(delta*rhosph_ep(ir))
1462            end do
1463 
1464 !          d2gam(1) = d2gam/drho+_drho+ (particle=ELECTRON), dgam/drho-_drho- (particle=POSITRON)
1465 !          d2gam(2) = d2gam/drho-_drho+ (particle=ELECTRON), dgam/drho+_drho- (particle=POSITRON)
1466 !          d2gam(3) = d2gam/drho-_drho- (particle=ELECTRON), dgam/drho+_drho+ (particle=POSITRON)
1467            ABI_ALLOCATE(d2gam,(mesh_size,2,3))
1468            d2gam(:,:,:)=zero
1469            do ir=1,mesh_size
1470              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
1471              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
1472              if (rho_ep_(ir)>tol14) then
1473                d2gam(ir,1,3)=(gam_(ir,1,4)+gam_(ir,1,5)-two*gam_(ir,1,1))/(delta*rho_ep_(ir))**2
1474                if (rho_(ir)>tol14) then
1475                  d2gam(ir,1,2)=(gam_(ir,1,6)+gam_(ir,1,7)+two*gam_(ir,1,1) &
1476 &                 -gam_(ir,1,2)-gam_(ir,1,3)-gam_(ir,1,4)-gam_(ir,1,5)) &
1477 &                 *half/(delta*rho_(ir))/(delta*rho_ep_(ir))
1478                end if
1479              end if
1480              if (rhosph_ep(ir)>tol14) then
1481                d2gam(ir,2,3)=(gam_(ir,2,4)+gam_(ir,2,5)-two*gam_(ir,2,1))/(delta*rhosph_ep(ir))**2
1482                if (rhosph(ir)>tol14) then
1483                  d2gam(ir,2,2)=(gam_(ir,2,6)+gam_(ir,2,7)+two*gam_(ir,2,1) &
1484 &                 -gam_(ir,2,2)-gam_(ir,2,3)-gam_(ir,2,4)-gam_(ir,2,5)) &
1485 &                 *half/(delta*rhosph(ir))/(delta*rhosph_ep(ir))
1486                end if
1487              end if
1488            end do
1489            ABI_DEALLOCATE(rho_)
1490            ABI_DEALLOCATE(rho_ep_)
1491 !          Compute useful sums of densities
1492            ABI_ALLOCATE(v1sum,(mesh_size,3))
1493            if ( dtset%pawxcdev>=2)  then
1494              ABI_ALLOCATE(v2sum,(mesh_size,lm_size,3))
1495            else
1496              ABI_ALLOCATE(v2sum,(0,0,0))
1497            end if
1498            rhotot(:,1)=sqfpi*rhosph(:);rhotot_ep(:,1)=sqfpi*rhosph_ep(:)
1499            call pawxcsum(1,1,1,lmselect,lmselect_ep,lm_size,mesh_size,3,dtset%pawxcdev,&
1500 &           pawang,rhotot,rhotot_ep,v1sum,v2sum)
1501 !          Compute final developpment of gamma moments
1502            gammam(:,:,:)=zero
1503            gammam(:,:,1)=gam_(:,:,1)*sqfpi
1504            gammam(:,1,1)=gammam(:,1,1)+(d2gam(:,1,2)*v1sum(:,2) &
1505 &           +half*(d2gam(:,1,1)*v1sum(:,1)+d2gam(:,1,3)*v1sum(:,3)))/sqfpi
1506            gammam(:,2,1)=gammam(:,2,1)+(d2gam(:,2,2)*v1sum(:,2) &
1507 &           +half*(d2gam(:,2,1)*v1sum(:,1)+d2gam(:,2,3)*v1sum(:,3)))/sqfpi
1508            do ilm=2,lm_size
1509              if (lmselect(ilm)) then
1510                gammam(:,1,ilm)=gammam(:,1,ilm)+d1gam(:,1,1)*rhotot(:,ilm)
1511                gammam(:,2,ilm)=gammam(:,2,ilm)+d1gam(:,2,1)*rhotot(:,ilm)
1512              end if
1513              if (lmselect_ep(ilm)) then
1514                gammam(:,1,ilm)=gammam(:,1,ilm)+d1gam(:,1,2)*rhotot_ep(:,ilm)
1515                gammam(:,2,ilm)=gammam(:,2,ilm)+d1gam(:,2,2)*rhotot_ep(:,ilm)
1516              end if
1517            end do
1518            if (dtset%pawxcdev>1) then
1519              do ilm=2,lm_size
1520                gammam(:,1,ilm)=gammam(:,1,ilm)+d2gam(:,1,2)*v2sum(:,ilm,2) &
1521 &               +half*(d2gam(:,1,1)*v2sum(:,ilm,1)+d2gam(:,1,3)*v2sum(:,ilm,3))
1522                gammam(:,2,ilm)=gammam(:,2,ilm)+d2gam(:,2,2)*v2sum(:,ilm,2) &
1523 &               +half*(d2gam(:,2,1)*v2sum(:,ilm,1)+d2gam(:,2,3)*v2sum(:,ilm,3))
1524              end do
1525            end if
1526            ABI_DEALLOCATE(gam_)
1527            ABI_DEALLOCATE(d1gam)
1528            ABI_DEALLOCATE(d2gam)
1529            ABI_DEALLOCATE(v1sum)
1530            ABI_DEALLOCATE(v2sum)
1531 
1532 !          Compute contribution to annihilation rate
1533 
1534 !          In state dependent scheme for Doppler, replace electronic density with one state
1535            if (option/=1) then
1536              rhotot(:,:) = rhotot_j(:,:)
1537              rhosph  (:) = rhosph_j  (:)
1538            end if
1539 
1540            ABI_ALLOCATE(gg,(mesh_size,4))
1541            gg=zero
1542            ABI_ALLOCATE(rhoarr1,(mesh_size))
1543            ABI_ALLOCATE(rhoarr2,(mesh_size))
1544            do ilm=1,lm_size
1545              do ilm1=1,lm_size
1546                if (lmselect(ilm1)) then
1547                  if (ilm1==1) rhoarr1(:)=sqfpi*rhosph(:)
1548                  if (ilm1/=1) rhoarr1(:)=rhotot(:,ilm1)
1549                  do ilm2=1,lm_size
1550                    if (lmselect_ep(ilm2)) then
1551                      if (ilm2==1) rhoarr2(:)=sqfpi*rhosph_ep(:)
1552                      if (ilm2/=1) rhoarr2(:)=rhotot_ep(:,ilm2)
1553                      if (ilm1>=ilm2) then
1554                        isel=pawang%gntselect(ilm,ilm2+ilm1*(ilm1-1)/2)
1555                      else
1556                        isel=pawang%gntselect(ilm,ilm1+ilm2*(ilm2-1)/2)
1557                      end if
1558                      if (isel>0) then
1559                        fact=pawang%realgnt(isel)
1560                        gg(:,1)=gg(:,1)+fact*rhoarr1(:)*rhoarr2(:)*gammam(:,1,ilm)
1561                        gg(:,2)=gg(:,2)+fact*rhoarr1(:)*rhoarr2(:)*gammam(:,2,ilm)
1562                      end if
1563                    end if
1564                  end do
1565                end if
1566              end do
1567            end do
1568            ABI_DEALLOCATE(rhoarr1)
1569            ABI_DEALLOCATE(rhoarr2)
1570            if (electronpositron%particle==EP_ELECTRON) then
1571              do ilm=1,lm_size
1572                if (lmselect(ilm)) gg(:,3)=gg(:,3)+rhotot(:,ilm)*rhocore(:)*gammam(:,1,ilm)
1573              end do
1574              gg(:,4)=sqfpi*rhotot(:,1)*rhocore(:)
1575            else if (electronpositron%particle==EP_POSITRON) then
1576              do ilm=1,lm_size
1577                if (lmselect_ep(ilm)) gg(:,3)=gg(:,3)+rhotot_ep(:,ilm)*rhocore(:)*gammam(:,1,ilm)
1578              end do
1579              gg(:,4)=sqfpi*rhotot_ep(:,1)*rhocore(:)
1580            end if
1581            do ii=1,4
1582              gg(1:mesh_size,ii)=gg(1:mesh_size,ii)*pawrad(itypat)%rad(1:mesh_size)**2
1583              call simp_gen(intg,gg(:,ii),pawrad(itypat))
1584              if (ii==1) lambda_paw         =lambda_paw         +lsign(iloop)*intg
1585              if (ii==2) lambda_paw_ipm     =lambda_paw_ipm     +lsign(iloop)*intg
1586              if (ii==3) lambda_core_paw    =lambda_core_paw    +lsign(iloop)*intg
1587              if (ii==4) lambda_core_paw_ipm=lambda_core_paw_ipm+lsign(iloop)*intg
1588            end do
1589            ABI_DEALLOCATE(gg)
1590            ABI_DEALLOCATE(gammam)
1591            ABI_DEALLOCATE(rhotot)
1592            ABI_DEALLOCATE(rhotot_ep)
1593            ABI_DEALLOCATE(rhosph)
1594            ABI_DEALLOCATE(rhosph_ep)
1595            if (option/=1) then
1596              ABI_DEALLOCATE(rhotot_j)
1597              ABI_DEALLOCATE(rhosph_j)
1598            end if
1599 
1600          end if ! dtset%pawxcdev
1601 
1602          ABI_DEALLOCATE(rhocore)
1603 
1604        end do ! iloop
1605 
1606        ABI_DEALLOCATE(rho1)
1607        ABI_DEALLOCATE(trho1)
1608        ABI_DEALLOCATE(rho1_ep)
1609        ABI_DEALLOCATE(trho1_ep)
1610        ABI_DEALLOCATE(rho1_j)
1611        ABI_DEALLOCATE(trho1_j)
1612        ABI_DEALLOCATE(nhat1)
1613        ABI_DEALLOCATE(nhat1_ep)
1614        ABI_DEALLOCATE(nhat1_j)
1615        ABI_DEALLOCATE(lmselect)
1616        ABI_DEALLOCATE(lmselect_ep)
1617        ABI_DEALLOCATE(lmselect_dum)
1618 
1619      end do ! iatom
1620 
1621 !    Reduction in case of distribution over atomic sites
1622      if (mpi_enreg%nproc_atom>1) then
1623        mpibuf(1)=lambda_paw     ;mpibuf(2)=lambda_paw_ipm
1624        mpibuf(3)=lambda_core_paw;mpibuf(4)=lambda_core_paw_ipm
1625        call xmpi_sum(mpibuf,mpi_enreg%comm_atom,ierr)
1626        lambda_paw=mpibuf(1)     ;lambda_paw_ipm=mpibuf(2)
1627        lambda_core_paw=mpibuf(3);lambda_core_paw_ipm=mpibuf(4)
1628      end if
1629 
1630 !    Add plane-wave and PAW contributions to annihilation rates
1631      lambda         =lambda         +lambda_paw
1632      lambda_ipm     =lambda_ipm     +lambda_paw_ipm
1633      lambda_core    =lambda_core    +lambda_core_paw
1634      lambda_core_ipm=lambda_core_ipm+lambda_core_paw_ipm
1635    end if ! dtset%usepaw
1636 
1637 
1638 !  Convert into proper units and print
1639 !  ---------------------------------------------------------------------------------------
1640 
1641 !  Sum valence and core contributions to annihilation rates
1642    lambda        =lambda        +lambda_core
1643    lambda_ipm    =lambda_ipm    +lambda_core_ipm
1644    if (dtset%usepaw==1) then
1645      lambda_paw    =lambda_paw    +lambda_core_paw
1646      lambda_paw_ipm=lambda_paw_ipm+lambda_core_paw_ipm
1647    end if
1648 
1649 !  Set annihilation rate in proper unit (picosec.)
1650    units=pi*(one/InvFineStruct)**3/Time_Sec/1.e12_dp/electronpositron%posocc
1651 
1652    lambda         =lambda         *units
1653    lambda_ipm     =lambda_ipm     *units
1654    lambda_core    =lambda_core    *units
1655    lambda_core_ipm=lambda_core_ipm*units
1656    lifetime    =one/lambda
1657    lifetime_ipm=one/lambda_ipm
1658    electronpositron%lambda=lambda
1659    electronpositron%lifetime=lifetime
1660    if (dtset%usepaw==1) then
1661      lambda_paw         =lambda_paw         *units
1662      lambda_paw_ipm     =lambda_paw_ipm     *units
1663      lambda_core_paw    =lambda_core_paw    *units
1664      lambda_core_paw_ipm=lambda_core_paw_ipm*units
1665    end if
1666    rate=lambda-lambda_core-lambda_paw+lambda_core_paw
1667    rate_paw=lambda_paw-lambda_core_paw
1668 !  Print life time and additional information
1669    if (option==1) then
1670      if (igam==1) then
1671        write(msg,'(a,80("-"),2a)') ch10,ch10,' Results for electron-positron annihilation:'
1672        call wrtout(ab_out,msg,'COLL')
1673        call wrtout(std_out,msg,'COLL')
1674      end if
1675      if (ngamma>1.and.igam==1) then
1676        write(msg,'(a,i1,a)') ch10,ngamma,&
1677 &       ' computations of positron lifetime have been performed (with different enhancement factors).'
1678        call wrtout(ab_out,msg,'COLL')
1679        call wrtout(std_out,msg,'COLL')
1680      end if
1681      if (ngamma>1) then
1682        write(msg,'(2a,i1)') ch10,"########## Lifetime computation ",igam
1683        call wrtout(ab_out,msg,'COLL')
1684        call wrtout(std_out,msg,'COLL')
1685      end if
1686      if (abs(electronpositron%ixcpositron)==1) then
1687        write(msg,'(4a)') ch10,' # Zero-positron density limit of Arponen and Pajanne provided by Boronski & Nieminen',&
1688 &       ch10,'   Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]]
1689      else if (electronpositron%ixcpositron==11) then
1690        write(msg,'(4a)') ch10,' # Zero-positron density limit of Arponen and Pajanne fitted by Sterne & Kaiser',&
1691 &       ch10,'   Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)' ! [[cite:Sterne1991]]
1692      else if (electronpositron%ixcpositron==2) then
1693        write(msg,'(4a)') ch10,' # Electron-positron correlation provided by Puska, Seitsonen, and Nieminen',&
1694 &       ch10,'   Ref: M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994)' !  [[cite:Puska1994]]
1695      else if (electronpositron%ixcpositron==3) then
1696        write(msg,'(8a)') ch10,' # Zero-positron density limit of Arponen and Pajanne provided by Boronski & Nieminen',&
1697 &       ch10,'   + GGA corrections',&
1698 &       ch10,'   Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)',& ! [[cite:Boronski1986]]
1699 &       ch10,'         B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1995)' ! [[cite:Barbiellini1995]]
1700      else if (electronpositron%ixcpositron==31) then
1701        write(msg,'(8a)') ch10,' # Zero-positron density limit of Arponen and Pajanne fitted by Sterne & Kaiser',&
1702 &       ch10,'   + GGA corrections',&
1703 &       ch10,'   Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)',& ! [[cite:Sterne1991]]
1704 &       ch10,'         B. Barbiellini, M.J. Puska, T. Torsti and R.M. Nieminen, Phys. Rev. B 51, 7341 (1995)' ! [[cite:Barbiellini1995]]
1705      end if
1706      call wrtout(ab_out,msg,'COLL')
1707      call wrtout(std_out,  msg,'COLL')
1708      if (igamma(igam)==0) then
1709        write(msg,'(a)')       ' # Enhancement factor set to one (test)'
1710      else if (igamma(igam)==1) then
1711        write(msg,'(3a)')      ' # Enhancement factor of Boronski & Nieminen',&
1712 &       ch10,'   Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]]
1713      else if (igamma(igam)==2) then
1714        write(msg,'(3a)')      ' # Enhancement factor of Boronski & Nieminen IN THE RPA LIMIT',&
1715 &       ch10,'   Ref.: Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986)' ! [[cite:Boronski1986]]
1716      else if (igamma(igam)==3) then
1717        write(msg,'(3a)')      ' # Enhancement factor of Sterne & Kaiser',&
1718 &       ch10,'   Ref.: P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991)' ! [[cite:Sterne1991]]
1719      else if (igamma(igam)==4) then
1720        write(msg,'(3a)')      ' # Enhancement factor of Puska, Seitsonen, and Nieminen',&
1721 &       ch10,'   Ref.: M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994)' !  [[cite:Puska1994]]
1722      end if
1723      call wrtout(ab_out,msg,'COLL')
1724      call wrtout(std_out,msg,'COLL')
1725      write(msg, '(4(2a,es16.8))' ) ch10,&
1726 &     ' Positron lifetime                         (ps)   =',lifetime    ,ch10,&
1727 &     ' Positron lifetime with IPM for core elec. (ps)   =',lifetime_ipm,ch10,&
1728 &     ' Annihilation rate                         (ns-1) =',lambda    *1000._dp,ch10,&
1729 &     ' Annihilation rate with IPM for core elec. (ns-1) =',lambda_ipm*1000._dp
1730      call wrtout(ab_out,msg,'COLL')
1731      call wrtout(std_out,msg,'COLL')
1732      write(msg,'(2a,5(2a,es16.8))' ) ch10,&
1733 &     ' Annihilation rate core/valence decomposition:',ch10,&
1734 &     '   Core    contribution to ann.rate          (ns-1) =', lambda_core                 *1000._dp,ch10,&
1735 &     '   Valence contribution to ann.rate          (ns-1) =',(lambda-lambda_core)         *1000._dp,ch10,&
1736 &     '   Core    contribution to ann.rate with IPM (ns-1) =', lambda_core_ipm             *1000._dp,ch10,&
1737 &     '   Valence contribution to ann.rate with IPM (ns-1) =',(lambda_ipm-lambda_core_ipm) *1000._dp
1738      call wrtout(ab_out,msg,'COLL')
1739      call wrtout(std_out,msg,'COLL')
1740      if (dtset%usepaw==1) then
1741        write(msg, '(2a,6(2a,es16.8))' ) ch10,&
1742 &       ' Annihilation rate PAW decomposition:',ch10,&
1743 &       '   Plane-wave contribution to ann.rate          (ns-1) =',(lambda-lambda_paw)*1000._dp,ch10,&
1744 &       '   Plane-wave valence contribution to ann.rate  (ns-1) =',(lambda-lambda_paw-lambda_core+lambda_core_paw)*1000._dp,ch10,&
1745 &       '   On-site core contribution to ann.rate        (ns-1) =', lambda_core_paw*1000._dp,ch10,&
1746 &       '   On-site valence contribution to ann.rate     (ns-1) =',(lambda_paw-lambda_core_paw)*1000._dp,ch10,&
1747 &       '   Plane-wave contribution to ann.rate with IPM (ns-1) =',(lambda_ipm-lambda_paw_ipm)*1000._dp,ch10,&
1748 &       '   Plane-wave core contrb. to ann.rate with IPM (ns-1) =',(lambda_core_ipm-lambda_core_paw_ipm)*1000._dp
1749        call wrtout(ab_out,msg,'COLL')
1750        call wrtout(std_out,msg,'COLL')
1751      end if
1752      if (dtset%usepaw==0.and.igam==ngamma) then ! These tests are not relevant with PAW
1753        write(msg, '(2a,3(2a,es16.8))' ) ch10,&
1754 &       ' ########## Some checks, for testing purpose:',ch10,&
1755 &       '   Number of core electrons      =',nbec,ch10,&
1756 &       '   Number of valence electrons   =',nbev,ch10,&
1757 &       '   Number of positrons           =',nbp
1758        call wrtout(ab_out,msg,'COLL')
1759        call wrtout(std_out,msg,'COLL')
1760      end if
1761    end if !end if option
1762  end do ! Big loop on igam
1763 
1764  if (option==1) then
1765    write(msg, '(3a)' ) ch10,'      (*) IPM=Independent particle Model',ch10
1766    call wrtout(ab_out,msg,'COLL')
1767    call wrtout(std_out,msg,'COLL')
1768  end if !end if option
1769 
1770 !Deallocate memory
1771  ABI_DEALLOCATE(igamma)
1772  if (dtset%usepaw==1.and.(.not.include_nhat_in_gamma)) then
1773    ABI_DEALLOCATE(rhor_)
1774    ABI_DEALLOCATE(rhor_ep_)
1775  end if
1776 
1777  DBG_EXIT("COLL")
1778 
1779 end subroutine poslifetime

ABINIT/posratecore [ Functions ]

[ Top ] [ 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

PARENTS

      posdoppler

CHILDREN

      gammapositron,mkdenpos,nderiv_gen,pawdensities,pawxcsum,simp_gen
      xmpi_sum

SOURCE

3358 subroutine posratecore(dtset,electronpositron,iatom,my_natom,mesh_sizej,mpi_enreg,&
3359 &                      option,pawang,pawrad,pawrhoij,pawrhoij_ep,&
3360 &                      pawtab,rate,rhocorej)
3361 
3362 
3363 !This section has been created automatically by the script Abilint (TD).
3364 !Do not modify the following lines by hand.
3365 #undef ABI_FUNC
3366 #define ABI_FUNC 'posratecore'
3367 !End of the abilint section
3368 
3369  implicit none
3370 
3371 !Arguments ------------------------------------
3372 !scalars
3373  integer,intent(in) :: iatom,my_natom,option,mesh_sizej
3374  real(dp),intent(out) :: rate
3375  type(dataset_type), intent(in) :: dtset
3376  type(electronpositron_type),pointer :: electronpositron
3377  type(MPI_type),intent(in) :: mpi_enreg
3378  type(pawang_type), intent(in) :: pawang
3379 !arrays
3380  real(dp),intent(in) :: rhocorej(mesh_sizej)
3381  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw)
3382  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*dtset%usepaw)
3383  type(pawrhoij_type),intent(in),target :: pawrhoij_ep(my_natom*dtset%usepaw)
3384  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
3385 
3386 !Local variables-------------------------------
3387 !scalars
3388  integer :: cplex,ierr,igamma,ii,ilm,ipt,ir
3389  integer :: itypat,iwarn,iwarnj,iwarnp,lm_size,lmn2_size,mesh_size
3390  integer :: ngr,ngrad,nspden_ep,opt_dens
3391  logical,parameter :: include_nhat_in_gamma=.false.
3392  real(dp),parameter :: delta=1.d-4
3393  real(dp) :: fact,fact2,intg
3394  real(dp) :: mpibuf,rdum,sqfpi
3395  character(len=500) :: msg
3396 !arrays
3397  logical,allocatable :: lmselect(:),lmselect_ep(:),lmselect_dum(:)
3398  real(dp),parameter :: qphon(3)=(/zero,zero,zero/),lsign(2)=(/one,-one/)
3399  real(dp),allocatable :: d1gam(:,:),d2gam(:,:),ff(:),gam_(:,:,:),gamma(:,:),gammam(:,:),gg(:,:)
3400  real(dp),allocatable :: grhocore2(:),grhocor2_(:),grhoe2(:),grho2_(:)
3401  real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:)
3402  real(dp),allocatable :: rho_(:),rho_ep_(:),rho1(:,:,:),rho1_ep(:,:,:)
3403  real(dp),allocatable :: rhoarr1(:),rhoarr1_ep(:),rhoarr2(:)
3404  real(dp),allocatable :: rhocore(:),rhocor_(:)
3405  real(dp),allocatable :: rhosph(:),rhosph_ep(:),rhotot(:,:),rhotot_ep(:,:)
3406  real(dp),allocatable :: trho1(:,:,:),trho1_ep(:,:,:)
3407  real(dp),allocatable :: v1sum(:,:),v2sum(:,:,:)
3408  type(pawrhoij_type),pointer :: pawrhoij_ep_(:)
3409 
3410 ! *************************************************************************
3411 
3412  DBG_ENTER("COLL")
3413 
3414 !Tests for developers
3415  if (.not.associated(electronpositron)) then
3416    msg='electronpositron variable must be associated!'
3417    MSG_BUG(msg)
3418  end if
3419 !Constants
3420  fact=0.0
3421  cplex=1;nspden_ep=1
3422  ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
3423  iwarn=0;iwarnj=0;iwarnp=1
3424  sqfpi=sqrt(four_pi)
3425 
3426 !Compatibility tests
3427  if (electronpositron%particle==EP_NOTHING) then
3428    msg='Not valid for electronpositron%particle=NOTHING!'
3429    MSG_BUG(msg)
3430  end if
3431 
3432  if (dtset%usepaw==1) then
3433    if(dtset%pawxcdev==0.and.ngrad==2) then
3434      msg='GGA is not implemented for pawxcdev=0 (use dtset%pawxcdev/=0)!'
3435      MSG_BUG(msg)
3436    end if
3437  end if
3438 
3439 !Select type(s) of enhancement factor
3440  if (electronpositron%ixcpositron==-1) igamma=0
3441  if (electronpositron%ixcpositron== 2) igamma=4
3442  if (electronpositron%ixcpositron==11.or.electronpositron%ixcpositron==31) igamma=3
3443  if (electronpositron%ixcpositron==1.or.electronpositron%ixcpositron==3) igamma=2
3444  if (option==2) igamma=0
3445 
3446  pawrhoij_ep_ => pawrhoij_ep
3447 
3448  rate=zero
3449 
3450  itypat=pawrhoij(iatom)%itypat
3451  lmn2_size=pawtab(itypat)%lmn2_size
3452  mesh_size=pawtab(itypat)%mesh_size
3453  lm_size=pawtab(itypat)%lcut_size**2
3454  cplex=1
3455  ngr=0;if (ngrad==2) ngr=mesh_size
3456 
3457 !Allocations of "on-site" densities
3458  ABI_ALLOCATE(rho1 ,(cplex*mesh_size,lm_size,nspden_ep))
3459  ABI_ALLOCATE(trho1,(cplex*mesh_size,lm_size,nspden_ep))
3460  ABI_ALLOCATE(rho1_ep ,(cplex*mesh_size,lm_size,nspden_ep))
3461  ABI_ALLOCATE(trho1_ep,(cplex*mesh_size,lm_size,nspden_ep))
3462  ABI_ALLOCATE(lmselect,(lm_size))
3463  ABI_ALLOCATE(lmselect_ep,(lm_size))
3464  ABI_ALLOCATE(lmselect_dum,(lm_size))
3465  if (include_nhat_in_gamma) then
3466    ABI_ALLOCATE(nhat1,(cplex*mesh_size,lm_size,nspden_ep))
3467    ABI_ALLOCATE(nhat1_ep,(cplex*mesh_size,lm_size,nspden_ep))
3468  else
3469    ABI_ALLOCATE(nhat1,(0,0,0))
3470    ABI_ALLOCATE(nhat1_ep,(0,0,0))
3471  end if
3472 
3473 !Compute "on-site" densities (n1, ntild1, nhat1) for electron and positron =====
3474  lmselect(:)=.true.
3475  opt_dens=1;if (include_nhat_in_gamma) opt_dens=0
3476  call pawdensities(rdum,cplex,iatom,lmselect,lmselect_dum,lm_size,nhat1,nspden_ep,1,&
3477 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij(iatom),&
3478 & pawtab(itypat),rho1,trho1)
3479  lmselect_ep(:)=.true.
3480  call pawdensities(rdum,cplex,iatom,lmselect_ep,lmselect_dum,lm_size,nhat1_ep,nspden_ep,1,&
3481 & 0,opt_dens,-1,0,pawang,0,pawrad(itypat),pawrhoij_ep_(iatom),&
3482 & pawtab(itypat),rho1_ep,trho1_ep)
3483 !Compute contribution to annihilation rate
3484 
3485  ABI_ALLOCATE(rhocore,(mesh_size))
3486 
3487 !First formalism: use densities on r,theta,phi
3488  if (dtset%pawxcdev==0) then
3489 
3490    ABI_ALLOCATE(gamma,(mesh_size,2))
3491    ABI_ALLOCATE(rhoarr1,(mesh_size))
3492    ABI_ALLOCATE(rhoarr1_ep,(mesh_size))
3493 
3494 !  Loop on the angular part
3495    do ipt=1,pawang%angl_size
3496 !    Build densities
3497      rhoarr1=zero;rhoarr1_ep=zero;rhocore=zero
3498      do ilm=1,lm_size
3499        if (lmselect(ilm)) rhoarr1(:)=rhoarr1(:)+rho1(:,ilm,1)*pawang%ylmr(ilm,ipt)
3500      end do
3501      do ilm=1,lm_size
3502        if (lmselect_ep(ilm)) rhoarr1_ep(:)=rhoarr1_ep(:)+rho1_ep(:,ilm,1)*pawang%ylmr(ilm,ipt)
3503      end do
3504      rhocore(:)=pawtab(itypat)%coredens(:)
3505 !    Make the densities positive
3506      if (electronpositron%particle==EP_ELECTRON) then
3507        call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1   ,dtset%xc_denpos)
3508        call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos)
3509      else if (electronpositron%particle==EP_POSITRON) then
3510        call mkdenpos(iwarn ,mesh_size,1,1,rhoarr1   ,dtset%xc_denpos)
3511        call mkdenpos(iwarnp,mesh_size,1,1,rhoarr1_ep,dtset%xc_denpos)
3512      end if
3513 !    Compute Gamma
3514      ABI_ALLOCATE(grhoe2,(ngr))
3515      ABI_ALLOCATE(grhocore2,(ngr))
3516      if (electronpositron%particle==EP_ELECTRON) then
3517        call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,mesh_size,&
3518 &       rhocore,rhoarr1_ep,rhoarr1,1)
3519      else if (electronpositron%particle==EP_POSITRON) then
3520        call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,mesh_size,&
3521 &       rhocore,rhoarr1,rhoarr1_ep,1)
3522      end if
3523      ABI_DEALLOCATE(grhoe2)
3524      ABI_DEALLOCATE(grhocore2)
3525 !    Compute contribution to annihilation rates
3526 
3527      ABI_ALLOCATE(ff,(mesh_size))
3528      ff(1:mesh_size)=rhoarr1_ep(1:mesh_size)*rhocorej(1:mesh_size) &
3529 &     *gamma(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
3530      call simp_gen(intg,ff,pawrad(itypat))
3531      intg=intg*pawang%angwgth(ipt)*four_pi
3532      rate         =rate         +intg
3533      ABI_DEALLOCATE(ff)
3534    end do ! ipt
3535    ABI_DEALLOCATE(gamma)
3536    ABI_DEALLOCATE(rhoarr1)
3537    ABI_DEALLOCATE(rhoarr1_ep)
3538 
3539 !Second formalism: use (l,m) moments for densities
3540  else if (dtset%pawxcdev/=0) then
3541 
3542 !  Build densities
3543    ABI_ALLOCATE(gammam,(mesh_size,lm_size))
3544    ABI_ALLOCATE(rhotot,(mesh_size,lm_size))
3545    ABI_ALLOCATE(rhotot_ep,(mesh_size,lm_size))
3546    ABI_ALLOCATE(rhosph,(mesh_size))
3547    ABI_ALLOCATE(rhosph_ep,(mesh_size))
3548 
3549    rhotot   (:,:)=rho1   (:,:,1)
3550    rhotot_ep(:,:)=rho1_ep(:,:,1)
3551    rhocore(:)=pawtab(itypat)%coredens(:)
3552    rhosph   (:)=rhotot   (:,1)/sqfpi
3553    rhosph_ep(:)=rhotot_ep(:,1)/sqfpi
3554 !  Make spherical densities positive
3555    if (electronpositron%particle==EP_ELECTRON) then
3556      call mkdenpos(iwarnp,mesh_size,1,1,rhosph   ,dtset%xc_denpos)
3557      call mkdenpos(iwarn ,mesh_size,1,1,rhosph_ep,dtset%xc_denpos)
3558    else if (electronpositron%particle==EP_POSITRON) then
3559      call mkdenpos(iwarn ,mesh_size,1,1,rhosph   ,dtset%xc_denpos)
3560      call mkdenpos(iwarnp,mesh_size,1,1,rhosph_ep,dtset%xc_denpos)
3561    end if
3562 
3563 !  Need gradients of electronic densities for GGA
3564    ABI_ALLOCATE(grhoe2,(ngr))
3565    ABI_ALLOCATE(grhocore2,(ngr))
3566    if (ngr>0) then
3567      if (electronpositron%particle==EP_ELECTRON) then
3568        call nderiv_gen(grhoe2,rhosph_ep,pawrad(itypat))
3569      else if (electronpositron%particle==EP_POSITRON) then
3570        call nderiv_gen(grhoe2,rhosph,pawrad(itypat))
3571      end if
3572      grhoe2(:)=grhoe2(:)**2
3573      call nderiv_gen(grhocore2,rhocore,pawrad(itypat))
3574      grhocore2(:)=grhocore2(:)**2
3575    end if
3576 !  Compute Gamma for (rho-,rho+),
3577 !  (rho- +drho-,rho+), (rho- -drho-,rho+),
3578 !  (rho-,rho+ +drho+), (rho-,rho+ -drho+),
3579 !  (rho- +drho-,rho+ +drho+), (rho- -drho-,rho+ -drho+)
3580 !  Do a seven steps loop
3581    ABI_ALLOCATE(gam_,(mesh_size,2,7))
3582    ABI_ALLOCATE(rho_,(mesh_size))
3583    ABI_ALLOCATE(rho_ep_,(mesh_size))
3584    ABI_ALLOCATE(rhocor_,(mesh_size))
3585    ABI_ALLOCATE(grho2_,(ngr))
3586    ABI_ALLOCATE(grhocor2_,(ngr))
3587 
3588    do ii=1,7
3589 !    Apply delta to get perturbed densities
3590      rho_(:)=rhosph(:);rho_ep_(:)=rhosph_ep(:);rhocor_(:)=rhocore(:)
3591      if (ngr>0) grho2_(:)=grhoe2(:)
3592      if (ngr>0) grhocor2_(:)=grhocore2(:)
3593      if (ii==2.or.ii==4.or.ii==6) fact=(one+delta)
3594      if (ii==3.or.ii==5.or.ii==7) fact=(one-delta)
3595      fact2=fact**2
3596      if (ii==2.or.ii==3.or.ii==6.or.ii==7) then
3597        rho_(:)=fact*rho_(:)
3598        if (electronpositron%particle==EP_POSITRON) then
3599          if (ngr>0) grho2_(:)=fact2*grho2_(:)
3600          rhocor_(:)=fact*rhocor_(:)
3601          if (ngr>0) grhocor2_(:)=fact2*grhocor2_(:)
3602        end if
3603      end if
3604 
3605      if (ii==4.or.ii==5.or.ii==6.or.ii==7) then
3606        rho_ep_(:)=fact*rho_ep_(:)
3607        if (electronpositron%particle==EP_ELECTRON) then
3608          if (ngr>0) grho2_(:)=fact2*grho2_(:)
3609          rhocor_(:)=fact*rhocor_(:)
3610          if (ngr>0) grhocor2_(:)=fact2*grhocor2_(:)
3611        end if
3612      end if
3613 !    Compute gamma for these perturbed densities
3614      if (electronpositron%particle==EP_ELECTRON) then
3615        call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma,ngr,mesh_size,rhocor_,rho_ep_,rho_,1)
3616      else if (electronpositron%particle==EP_POSITRON) then
3617        call gammapositron(gam_(:,:,ii),grhocor2_,grho2_,igamma,ngr,mesh_size,rhocor_,rho_,rho_ep_,1)
3618      end if
3619 
3620    end do ! end loop ii=1,7
3621 
3622    ABI_DEALLOCATE(rhocor_)
3623    ABI_DEALLOCATE(grho2_)
3624    ABI_DEALLOCATE(grhocor2_)
3625    ABI_DEALLOCATE(grhoe2)
3626    ABI_DEALLOCATE(grhocore2)
3627    rho_   (:)=rhosph   (:);if (electronpositron%particle==EP_POSITRON) rho_   (:)=rho_   (:)+rhocore(:)
3628    rho_ep_(:)=rhosph_ep(:);if (electronpositron%particle==EP_ELECTRON) rho_ep_(:)=rho_ep_(:)+rhocore(:)
3629 !  Compute numerical first and second derivatives of Gamma
3630 !  d1gam(1) = dgam/drho+ (particle=ELECTRON), dgam/drho- (particle=POSITRON)
3631 !  d1gam(2) = dgam/drho- (particle=ELECTRON), dgam/drho+ (particle=POSITRON)
3632    ABI_ALLOCATE(d1gam,(mesh_size,2))
3633    d1gam(:,:)=zero
3634    do ir=1,mesh_size
3635      if (rho_     (ir)>tol14) d1gam(ir,1)=(gam_(ir,1,2)-gam_(ir,1,3))*half/(delta*rho_     (ir))
3636      if (rho_ep_  (ir)>tol14) d1gam(ir,2)=(gam_(ir,1,4)-gam_(ir,1,5))*half/(delta*rho_ep_  (ir))
3637    end do
3638 
3639 !  d2gam(1) = d2gam/drho+_drho+ (particle=ELECTRON), dgam/drho-_drho- (particle=POSITRON)
3640 !  d2gam(2) = d2gam/drho-_drho+ (particle=ELECTRON), dgam/drho+_drho- (particle=POSITRON)
3641 !  d2gam(3) = d2gam/drho-_drho- (particle=ELECTRON), dgam/drho+_drho+ (particle=POSITRON)
3642    ABI_ALLOCATE(d2gam,(mesh_size,3))
3643    d2gam(:,:)=zero
3644    do ir=1,mesh_size
3645      if (rho_  (ir)>tol14) d2gam(ir,1)=(gam_(ir,1,2)+gam_(ir,1,3)-two*gam_(ir,1,1))/(delta*rho_  (ir))**2
3646      if (rho_ep_(ir)>tol14) then
3647        d2gam(ir,3)=(gam_(ir,1,4)+gam_(ir,1,5)-two*gam_(ir,1,1))/(delta*rho_ep_(ir))**2
3648        if (rho_(ir)>tol14) then
3649          d2gam(ir,2)=(gam_(ir,1,6)+gam_(ir,1,7)+two*gam_(ir,1,1) &
3650 &         -gam_(ir,1,2)-gam_(ir,1,3)-gam_(ir,1,4)-gam_(ir,1,5)) &
3651 &         *half/(delta*rho_(ir))/(delta*rho_ep_(ir))
3652        end if
3653      end if
3654    end do
3655 
3656    ABI_DEALLOCATE(rho_)
3657    ABI_DEALLOCATE(rho_ep_)
3658 !  Compute useful sums of densities
3659    ABI_ALLOCATE(v1sum,(mesh_size,3))
3660    if ( dtset%pawxcdev>=2)  then
3661      ABI_ALLOCATE(v2sum,(mesh_size,lm_size,3))
3662    else
3663      ABI_ALLOCATE(v2sum,(0,0,0))
3664    end if
3665    rhotot(:,1)=sqfpi*rhosph(:);rhotot_ep(:,1)=sqfpi*rhosph_ep(:)
3666    call pawxcsum(1,1,1,lmselect,lmselect_ep,lm_size,mesh_size,3,dtset%pawxcdev,&
3667 &   pawang,rhotot,rhotot_ep,v1sum,v2sum)
3668 !  Compute final developpment of gamma moments
3669    gammam(:,:)=zero
3670    gammam(:,1)=gam_(:,1,1)*sqfpi
3671    gammam(:,1)=gammam(:,1)+(d2gam(:,2)*v1sum(:,2) &
3672 &   +half*(d2gam(:,1)*v1sum(:,1)+d2gam(:,3)*v1sum(:,3)))/sqfpi
3673    do ilm=2,lm_size
3674      if (lmselect(ilm)) then
3675        gammam(:,ilm)=gammam(:,ilm)+d1gam(:,1)*rhotot(:,ilm)
3676      end if
3677      if (lmselect_ep(ilm)) then
3678        gammam(:,ilm)=gammam(:,ilm)+d1gam(:,2)*rhotot_ep(:,ilm)
3679      end if
3680    end do
3681    if (dtset%pawxcdev>1) then
3682      do ilm=2,lm_size
3683        gammam(:,ilm)=gammam(:,ilm)+d2gam(:,2)*v2sum(:,ilm,2) &
3684 &       +half*(d2gam(:,1)*v2sum(:,ilm,1)+d2gam(:,3)*v2sum(:,ilm,3))
3685      end do
3686    end if
3687 
3688    ABI_DEALLOCATE(gam_)
3689    ABI_DEALLOCATE(d1gam)
3690    ABI_DEALLOCATE(d2gam)
3691    ABI_DEALLOCATE(v1sum)
3692    ABI_DEALLOCATE(v2sum)
3693 !  Compute contribution to annihilation rate
3694    ABI_ALLOCATE(gg,(mesh_size,4))
3695    gg=zero
3696    ABI_ALLOCATE(rhoarr2,(mesh_size))
3697    do ilm=1,lm_size
3698      if (lmselect_ep(ilm)) gg(:,1)=gg(:,1)+rhotot_ep(:,ilm)*rhocorej(:)*gammam(:,ilm)
3699    end do
3700    ABI_DEALLOCATE(rhoarr2)
3701    gg(1:mesh_size,1)=gg(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
3702    call simp_gen(intg,gg(:,1),pawrad(itypat))
3703    rate         =rate         +intg
3704    ABI_DEALLOCATE(gg)
3705    ABI_DEALLOCATE(gammam)
3706    ABI_DEALLOCATE(rhotot)
3707    ABI_DEALLOCATE(rhotot_ep)
3708    ABI_DEALLOCATE(rhosph)
3709    ABI_DEALLOCATE(rhosph_ep)
3710 
3711  end if ! dtset%pawxcdev
3712  ABI_DEALLOCATE(rhocore)
3713 
3714  ABI_DEALLOCATE(rho1)
3715  ABI_DEALLOCATE(trho1)
3716  ABI_DEALLOCATE(rho1_ep)
3717  ABI_DEALLOCATE(trho1_ep)
3718  ABI_DEALLOCATE(lmselect)
3719  ABI_DEALLOCATE(lmselect_ep)
3720  ABI_DEALLOCATE(lmselect_dum)
3721  ABI_DEALLOCATE(nhat1)
3722  ABI_DEALLOCATE(nhat1_ep)
3723 
3724 !Reduction in case of distribution over atomic sites
3725  if (mpi_enreg%nproc_atom>1) then
3726    mpibuf=rate
3727    call xmpi_sum(mpibuf,mpi_enreg%comm_atom,ierr)
3728    rate=mpibuf
3729  end if
3730 
3731  DBG_EXIT("COLL")
3732 
3733 end subroutine posratecore

ABINIT/setup_positron [ Functions ]

[ Top ] [ 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
  fred(3,natom)=forces 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)
  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=informations 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
  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
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
  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)

PARENTS

      scfcv

CHILDREN

      energies_copy,energies_init,forstr,fourdp,getcut,hartre,initrhoij
      initro,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawmknhat,pawrhoij_alloc
      pawrhoij_copy,pawrhoij_free,read_rhor,wrtout

SOURCE

184 subroutine setup_positron(atindx,atindx1,cg,cprj,dtefield,dtfil,dtset,ecore,eigen,etotal,electronpositron,&
185 &          energies,fock,forces_needed,fred,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,hdr,ifirst_gs,indsym,istep,istep_mix,kg,&
186 &          kxc,maxfor,mcg,mcprj,mgfft,mpi_enreg,my_natom,n3xccc,nattyp,nfft,ngfft,ngrvdw,nhat,nkxc,npwarr,nvresid,occ,optres,&
187 &          paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1dc,psps,rhog,rhor,&
188 &          rprimd,stress_needed,strsxc,symrec,ucvol,usecprj,vhartr,vpsp,vxc,&
189 &          xccc3d,xred,ylm,ylmgr)
190 
191 
192 !This section has been created automatically by the script Abilint (TD).
193 !Do not modify the following lines by hand.
194 #undef ABI_FUNC
195 #define ABI_FUNC 'setup_positron'
196 !End of the abilint section
197 
198  implicit none
199 
200 !Arguments ------------------------------------
201 !scalars
202  integer,intent(in) :: forces_needed,ifirst_gs,istep,mcg,mcprj,mgfft,my_natom,n3xccc,nfft
203  integer,intent(in) :: ngrvdw,nkxc,optres,stress_needed,usecprj
204  integer,intent(inout) :: istep_mix
205  real(dp),intent(in) :: ecore,etotal,gsqcut,maxfor,ucvol
206  type(efield_type),intent(in) :: dtefield
207  type(datafiles_type),intent(in) :: dtfil
208  type(dataset_type),intent(in) :: dtset
209  type(electronpositron_type),pointer :: electronpositron
210  type(energies_type),intent(inout) :: energies
211  type(hdr_type),intent(inout) :: hdr
212  type(MPI_type),intent(inout) :: mpi_enreg
213  type(pawang_type),intent(in) :: pawang
214  type(pawfgr_type),intent(in) :: pawfgr
215  type(pseudopotential_type), intent(in) :: psps
216 type(fock_type),pointer, intent(inout) :: fock
217 !arrays
218  integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
219  integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%natom),ngfft(18)
220  integer,intent(in) :: npwarr(dtset%nkpt),symrec(3,3,dtset%nsym)
221  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),grchempottn(3,dtset%natom)
222  real(dp),intent(in) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc)
223  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),ph1dc(2,(3*(2*dtset%mgfft+1)*dtset%natom)*dtset%usepaw)
224  real(dp),intent(in) :: rprimd(3,3),strsxc(6),vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden)
225  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
226  real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
227  real(dp),intent(inout) :: cg(2,mcg)
228  real(dp),intent(inout) :: nhat(nfft,dtset%nspden*dtset%usepaw)
229  real(dp),intent(inout) :: nvresid(nfft,dtset%nspden)
230  real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol),fred(3,dtset%natom)
231  real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
232  real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,dtset%nspden)
233  real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom)
234  type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj)
235  type(paw_ij_type),intent(in) :: paw_ij(my_natom*dtset%usepaw)
236  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*dtset%usepaw)
237  type(pawrad_type),intent(in)  :: pawrad(dtset%ntypat*dtset%usepaw)
238  type(pawtab_type),intent(in)  :: pawtab(dtset%ntypat*dtset%usepaw)
239  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*dtset%usepaw)
240 
241 !Local variables-------------------------------
242 !scalars
243  integer,parameter :: cplex1=1
244  integer :: history_level,iatom,iband,icalctype,icalctype0,icg,ifft,ikpt
245  integer :: iocc,ireadwf,ispden,isppol,n3xccc0,nocc,optfor,optstr,rdwrpaw,comm_cell
246  logical,parameter :: always_restart=.false.  ! Set to true to restart by a pure electronic step at each new atomic structure
247  logical :: need_scocc,new_calctype
248  real(dp) :: boxcut_dum,diffor_dum,ecut_eff,eigtmp,etotal_read,gsqcut_eff,maxfor_dum
249  real(dp) :: maxocc,nelect,occlast,occtmp,rhotmp
250  character(len=69) :: TypeCalcStrg
251  character(len=500) :: message
252  character(len=fnlen) :: fname
253  type(energies_type) :: energies_tmp
254  type(wvl_data) :: wvl
255  type(hdr_type) :: hdr_den
256 !arrays
257  integer,allocatable :: nlmn(:)
258  real(dp) :: cgtmp(2)
259  real(dp),parameter :: qphon(3)=(/zero,zero,zero/)
260  real(dp),allocatable :: favg_dum(:),fcart_dum(:,:),forold_dum(:,:),fred_tmp(:,:)
261  real(dp),allocatable :: gresid_dum(:,:),grhf_dum(:,:),grxc_dum(:,:)
262  real(dp),allocatable :: rhog_ep(:,:),scocc(:),str_tmp(:),synlgr_dum(:,:)
263  real(dp) :: nhatgr(0,0,0)
264  type(pawcprj_type),allocatable :: cprj_tmp(:,:)
265  type(pawrhoij_type),allocatable :: pawrhoij_tmp(:)
266 
267 ! *************************************************************************
268 
269 !Compatibility tests
270  if (dtset%positron==0) then
271    MSG_BUG('Not valid for dtset%positron=0!')
272  end if
273 
274  if (istep>1.and.nfft/=electronpositron%nfft) then
275    MSG_BUG('Invalid value for nfft!')
276  end if
277 
278  if (dtset%usewvl==1) then
279    MSG_BUG('Not valid for wavelets!')
280  end if
281 
282  if (dtset%positron==1) then
283    do isppol=1,dtset%nsppol
284      do ikpt=1,dtset%nkpt
285        if (dtset%nband(ikpt+dtset%nkpt*(isppol-1))/=dtset%nband(1)) then
286          message = "dtset%positron needs nband to be the same at each k-point !"
287          MSG_ERROR(message)
288        end if
289      end do
290    end do
291  end if
292 
293  comm_cell = mpi_enreg%comm_cell
294 
295 !-----------------------------------------------------------------------
296 !Compute new value for calctype (type of electron-positron calculation)
297 !-----------------------------------------------------------------------
298  icalctype0=electronpositron%calctype
299 
300  new_calctype=.false.
301  if (dtset%positron==1.or.dtset%positron==2) then
302    if (istep==1) new_calctype=.true.
303  else if (dtset%positron<0) then
304    if (ifirst_gs/=0.and.istep==1.and.(.not.always_restart)) new_calctype=.true.
305    if (electronpositron%scf_converged) new_calctype=.true.
306  end if
307 
308 !Comment:
309 !history_level=-1:  not used
310 !history_level= 0:  rhor from scratch, rhor_ep from scratch or read
311 !history_level= 1:  rhor in memory, rhor_ep from scratch or read
312 !history_level= 2:  rhor_ep <-rhor, rhor from scratch
313 !history_level= 3:  rhor_ep <-> rhor
314 !history_level= 4:  rhor in memory, rhor_ep in memory
315  history_level=-1
316  if (dtset%positron==1.or.dtset%positron==2) then
317    if (ifirst_gs==0.and.istep==1) history_level=0
318    if (ifirst_gs/=0.and.istep==1) history_level=4
319  else if (dtset%positron<0) then
320    if (.not.electronpositron%scf_converged) then
321      if (ifirst_gs/=0.and.istep==1.and.(.not.always_restart)) history_level=4
322    else if (electronpositron%scf_converged) then
323      if (icalctype0==0) history_level=2
324      if (icalctype0> 0) history_level=3
325    end if
326  end if
327 
328  electronpositron%calctype=icalctype0
329  if (dtset%positron==1.or.dtset%positron==2) then
330    electronpositron%calctype=dtset%positron
331  else if (dtset%positron<0) then
332    if (electronpositron%scf_converged) then
333      if (icalctype0==0) electronpositron%calctype=1
334      if (icalctype0>0 ) electronpositron%calctype=3-electronpositron%calctype
335    else if (ifirst_gs/=0.and.istep==1) then
336      if (always_restart) then
337        electronpositron%calctype=0
338      else
339        electronpositron%calctype=2
340 !       if (electronpositron%particle==EP_POSITRON) electronpositron%calctype=1
341      end if
342    end if
343  end if
344 
345  electronpositron%scf_converged=.false.
346  if (new_calctype) electronpositron%istep=electronpositron%istep+1
347  if (istep==1) electronpositron%istep=1
348  ireadwf=dtfil%ireadwf;if (electronpositron%istep>1) ireadwf=0
349 
350 !============================================
351 !The following lines occur only when the type
352 !of electron-positron calculation changes
353 !============================================
354  if (new_calctype) then
355 
356 !  Reset some indexes
357    if (electronpositron%calctype==0) then
358      electronpositron%particle=EP_NOTHING
359    else if (electronpositron%calctype==1) then
360      electronpositron%particle=EP_ELECTRON
361    else if (electronpositron%calctype==2) then
362      electronpositron%particle=EP_POSITRON
363    end if
364    electronpositron%has_pos_ham=mod(electronpositron%calctype,2)
365    electronpositron%istep_scf=1
366    istep_mix=1
367 
368 !  -----------------------------------------------------------------------------------------
369 !  Update forces and stresses
370 !  If electronpositron%calctype==1: fred_ep/stress_ep are the electronic fred/stress
371 !  If electronpositron%calctype==2: fred_ep/stress_ep are the positronic fred/stress
372 !  -----------------------------------------------------------------------------------------
373    if (history_level==2.or.history_level==3) then
374      optstr=0;optfor=0
375      if (allocated(electronpositron%stress_ep)) optstr=stress_needed
376      if (allocated(electronpositron%fred_ep).and.forces_needed==2) optfor=1
377      if (optfor>0.or.optstr>0) then
378        ABI_ALLOCATE(favg_dum,(3))
379        ABI_ALLOCATE(fcart_dum,(3,dtset%natom))
380        ABI_ALLOCATE(forold_dum,(3,dtset%natom))
381        ABI_ALLOCATE(gresid_dum,(3,dtset%natom))
382        ABI_ALLOCATE(grhf_dum,(3,dtset%natom))
383        ABI_ALLOCATE(grxc_dum,(3,dtset%natom))
384        ABI_ALLOCATE(synlgr_dum,(3,dtset%natom))
385        ABI_ALLOCATE(fred_tmp,(3,dtset%natom))
386        ABI_ALLOCATE(str_tmp,(6))
387        forold_dum=zero;n3xccc0=n3xccc
388        icalctype=electronpositron%calctype;electronpositron%calctype=-icalctype0
389        if (electronpositron%calctype==0) electronpositron%calctype=-100
390        if (electronpositron%calctype==-1) n3xccc0=0  ! Note: if calctype=-1, previous calculation was positron
391        call forstr(atindx1,cg,cprj,diffor_dum,dtefield,dtset,eigen,electronpositron,energies,&
392 &       favg_dum,fcart_dum,fock,forold_dum,fred_tmp,grchempottn,gresid_dum,grewtn,grhf_dum,grvdw,grxc_dum,gsqcut,&
393 &       indsym,kg,kxc,maxfor_dum,mcg,mcprj,mgfft,mpi_enreg,my_natom,n3xccc0,nattyp,nfft,ngfft,&
394 &       ngrvdw,nhat,nkxc,npwarr,dtset%ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,&
395 &       pawfgrtab,pawrad,pawrhoij,pawtab,ph1dc,ph1d,psps,rhog,rhor,rprimd,optstr,strsxc,str_tmp,symrec,&
396 &       synlgr_dum,ucvol,usecprj,vhartr,vpsp,vxc,wvl,xccc3d,xred,ylm,ylmgr,0.0_dp)
397        electronpositron%calctype=icalctype
398        if (optfor>0) electronpositron%fred_ep(:,:)=fred_tmp(:,:)
399        if (optstr>0) electronpositron%stress_ep(:)=str_tmp(:)
400        ABI_DEALLOCATE(favg_dum)
401        ABI_DEALLOCATE(fcart_dum)
402        ABI_DEALLOCATE(forold_dum)
403        ABI_DEALLOCATE(gresid_dum)
404        ABI_DEALLOCATE(grhf_dum)
405        ABI_DEALLOCATE(grxc_dum)
406        ABI_DEALLOCATE(synlgr_dum)
407        ABI_DEALLOCATE(fred_tmp)
408        ABI_DEALLOCATE(str_tmp)
409      end if
410      if (optfor==0.and.forces_needed>0.and.allocated(electronpositron%fred_ep)) then
411        electronpositron%fred_ep(:,:)=fred(:,:)-electronpositron%fred_ep(:,:)
412      end if
413    end if
414 
415 !  ----------------------------------------------------------------------------------------------------
416 !  Initialize/Update densities
417 !  If electronpositron%calctype==1: rhor is the positronic density, rhor_ep is the electronic density
418 !  If electronpositron%calctype==2: rhor is the electronic density, rhor_ep is the positronic density
419 !  ---------------------------------------------------------------------------------------------------
420    ABI_ALLOCATE(rhog_ep,(2,nfft))
421 
422 !  ===== PREVIOUS DENSITY RHOR_EP:
423    if (history_level==0.or.history_level==1) then
424 !    ----- Read from disk
425      if (dtset%positron>0) then
426        rdwrpaw=dtset%usepaw
427        fname=trim(dtfil%fildensin);if (dtset%positron==2) fname=trim(dtfil%fildensin)//'_POSITRON'
428        call read_rhor(trim(fname), cplex1, dtset%nspden, nfft, ngfft, rdwrpaw, mpi_enreg, electronpositron%rhor_ep, &
429        hdr_den, electronpositron%pawrhoij_ep, comm_cell, check_hdr=hdr)
430        etotal_read = hdr_den%etot; call hdr_free(hdr_den)
431        call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
432        if (dtset%usepaw==1.and.allocated(electronpositron%nhat_ep)) then
433          call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,&
434 &         dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,electronpositron%nhat_ep,&
435 &         electronpositron%pawrhoij_ep,electronpositron%pawrhoij_ep,pawtab,&
436 &         qphon,rprimd,ucvol,dtset%usewvl,xred,distribfft=mpi_enreg%distribfft,&
437 &         comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
438 &         comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0)
439        end if
440      end if
441 !    ----- Electronic from scratch
442      if (dtset%positron<0.and.electronpositron%calctype==1) then
443        ecut_eff=dtset%pawecutdg*(dtset%dilatmx)**2
444        call getcut(boxcut_dum,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,std_out,qphon,ngfft)
445        call initro(atindx,dtset%densty,gmet,gsqcut_eff,dtset%usepaw,mgfft,mpi_enreg,&
446 &       psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,dtset%ntypat,&
447 &       dtset%paral_kgb,psps,pawtab,ph1d,psps%qgrid_vl,rhog_ep,electronpositron%rhor_ep,&
448 &       dtset%spinat,ucvol,dtset%usepaw,dtset%ziontypat,dtset%znucl)
449        if (dtset%usepaw==1) then
450          if (size(electronpositron%pawrhoij_ep)>0) then
451            ABI_DATATYPE_ALLOCATE(pawrhoij_tmp,(my_natom))
452            call initrhoij(electronpositron%pawrhoij_ep(1)%cplex,dtset%lexexch,&
453 &           dtset%lpawu,my_natom,dtset%natom,dtset%nspden,&
454 &           electronpositron%pawrhoij_ep(1)%nspinor,dtset%nsppol,&
455 &           dtset%ntypat,pawrhoij_tmp,dtset%pawspnorb,pawtab,dtset%spinat,dtset%typat,&
456 &           ngrhoij=electronpositron%pawrhoij_ep(1)%ngrhoij,&
457 &           nlmnmix=electronpositron%pawrhoij_ep(1)%lmnmix_sz,&
458 &           use_rhoij_=electronpositron%pawrhoij_ep(1)%use_rhoij_,&
459 &           use_rhoijres=electronpositron%pawrhoij_ep(1)%use_rhoijres,&
460 &           comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
461            if (electronpositron%pawrhoij_ep(1)%lmnmix_sz>0) then
462              do iatom=1,my_natom
463                pawrhoij_tmp(iatom)%kpawmix(:)=electronpositron%pawrhoij_ep(iatom)%kpawmix(:)
464              end do
465            end if
466            call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep)
467            call pawrhoij_free(pawrhoij_tmp)
468            ABI_DATATYPE_DEALLOCATE(pawrhoij_tmp)
469          end if
470          if (allocated(electronpositron%nhat_ep)) then
471            call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,&
472 &           dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,electronpositron%nhat_ep,&
473 &           electronpositron%pawrhoij_ep,electronpositron%pawrhoij_ep,pawtab,qphon,rprimd,&
474 &           ucvol,dtset%usewvl,xred,distribfft=mpi_enreg%distribfft,&
475 &           comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
476 &           comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0)
477          end if
478        end if
479      end if
480 !    ----- Positronic from scratch
481      if (dtset%positron<0.and.electronpositron%calctype==2) then
482        electronpositron%rhor_ep(:,1)=one/ucvol
483        if (dtset%nspden>=2) electronpositron%rhor_ep(:,2)=half/ucvol
484        if (dtset%nspden==4) electronpositron%rhor_ep(:,3:4)=zero
485        rhog_ep=zero;rhog_ep(1,1)=one/ucvol
486        if (dtset%usepaw==1) then
487          do iatom=1,dtset%natom
488            electronpositron%pawrhoij_ep(iatom)%rhoijp(:,:)=zero
489            electronpositron%pawrhoij_ep(iatom)%nrhoijsel=0
490          end do
491          if (allocated(electronpositron%nhat_ep)) electronpositron%nhat_ep=zero
492        end if
493      end if
494    end if
495 !  ----- Deduced from rhor in memory
496    if (history_level==2) then
497      electronpositron%rhor_ep(:,:)=rhor(:,:)
498      rhog_ep(:,:)=rhog(:,:)
499      if (dtset%usepaw==1) then
500        call pawrhoij_copy(pawrhoij,electronpositron%pawrhoij_ep)
501        if (allocated(electronpositron%nhat_ep)) electronpositron%nhat_ep(:,:)=nhat(:,:)
502      end if
503    end if
504 
505 !  ===== CURRENT DENSITY RHOR:
506    if (history_level==0.or.history_level==2) then
507      if (ireadwf==0) then
508 !      ----- Positronic from scratch
509        if (electronpositron%calctype==1) then
510          rhor(:,1)=one/ucvol
511          if (dtset%nspden>=2) rhor(:,2)=half/ucvol
512          if (dtset%nspden==4) rhor(:,3:4)=zero
513          rhog=zero;rhog(1,1)=one/ucvol
514          if (dtset%usepaw==1) then
515            do iatom=1,my_natom
516              pawrhoij(iatom)%rhoijp(:,:)=zero
517              pawrhoij(iatom)%nrhoijsel=0
518            end do
519            nhat(:,:)=zero
520          end if
521        end if
522 !      ----- Electronic from scratch
523        if (electronpositron%calctype==2) then
524          ecut_eff=dtset%pawecutdg*(dtset%dilatmx)**2
525          call getcut(boxcut_dum,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,std_out,qphon,ngfft)
526          call initro(atindx,dtset%densty,gmet,gsqcut_eff,dtset%usepaw,mgfft,mpi_enreg,&
527 &         psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,dtset%ntypat,&
528 &         dtset%paral_kgb,psps,pawtab,ph1d,psps%qgrid_vl,rhog,rhor,dtset%spinat,ucvol,&
529 &         dtset%usepaw,dtset%ziontypat,dtset%znucl)
530 
531          if (dtset%usepaw==1) then
532            if (size(pawrhoij)>0) then
533              ABI_DATATYPE_ALLOCATE(pawrhoij_tmp,(my_natom))
534              call initrhoij(pawrhoij(1)%cplex,dtset%lexexch,dtset%lpawu,&
535 &             my_natom,dtset%natom,dtset%nspden,pawrhoij(1)%nspinor,dtset%nsppol,&
536 &             dtset%ntypat,pawrhoij_tmp,dtset%pawspnorb,pawtab,dtset%spinat,&
537 &             dtset%typat,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,&
538 &             use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres,&
539 &             comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
540              do iatom=1,my_natom
541                pawrhoij_tmp(iatom)%kpawmix(:)=pawrhoij(iatom)%kpawmix(:)
542              end do
543              call pawrhoij_copy(pawrhoij_tmp,pawrhoij)
544              call pawrhoij_free(pawrhoij_tmp)
545              ABI_DATATYPE_DEALLOCATE(pawrhoij_tmp)
546            end if
547            call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,&
548 &           dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,nhat,&
549 &           pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, &
550 &           comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
551 &           comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,&
552 &           me_g0=mpi_enreg%me_g0,distribfft=mpi_enreg%distribfft)
553          end if
554 
555        end if
556      end if
557    end if
558 
559 !  ===== EXCHANGE POSITRONIC AND ELECTRONIC DENSITY (CURRENT AND PREVIOUS)
560    if (history_level==3) then
561      do ispden=1,dtset%nspden
562        do ifft=1,nfft
563          rhotmp=rhor(ifft,ispden)
564          rhor(ifft,ispden)=electronpositron%rhor_ep(ifft,ispden)
565          electronpositron%rhor_ep(ifft,ispden)=rhotmp
566        end do
567      end do
568      rhog_ep(:,:)=rhog
569      call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
570 !    If PAW, exchange "positronic" and "electronic" rhoij
571      if (dtset%usepaw==1) then
572        if (size(pawrhoij)>0.and.size(electronpositron%pawrhoij_ep)>0) then
573          ABI_DATATYPE_ALLOCATE(pawrhoij_tmp,(my_natom))
574          call pawrhoij_alloc(pawrhoij_tmp,pawrhoij(1)%cplex,pawrhoij(1)%nspden,&
575 &         pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,dtset%typat,&
576 &         pawtab=pawtab,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,&
577 &         use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres, &
578 &         comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
579          call pawrhoij_copy(pawrhoij,pawrhoij_tmp)
580          call pawrhoij_copy(electronpositron%pawrhoij_ep,pawrhoij)
581          call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep)
582          call pawrhoij_free(pawrhoij_tmp)
583          ABI_DATATYPE_DEALLOCATE(pawrhoij_tmp)
584        end if
585        if (allocated(electronpositron%nhat_ep)) then
586          do ispden=1,dtset%nspden
587            do ifft=1,nfft
588              rhotmp=nhat(ifft,ispden)
589              nhat(ifft,ispden)=electronpositron%nhat_ep(ifft,ispden)
590              electronpositron%nhat_ep(ifft,ispden)=rhotmp
591            end do
592          end do
593        end if
594      end if
595    end if
596 
597 !  ===== COMPUTE HARTREE POTENTIAL ASSOCIATED TO RHOR_EP
598    if (history_level==4) then
599      call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
600    end if
601    if (history_level/=-1) then
602      call hartre(1,gsqcut,dtset%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog_ep,rprimd,&
603 &     electronpositron%vha_ep)
604      electronpositron%vha_ep=-electronpositron%vha_ep
605    else
606      electronpositron%vha_ep=zero
607    end if
608    ABI_DEALLOCATE(rhog_ep)
609 
610 !  ----------------------------------------------------------------------
611 !  Initialize/Update energies
612 !  ----------------------------------------------------------------------
613    electronpositron%etotal_prev=etotal
614    electronpositron%maxfor_prev=maxfor
615 
616 !  Inits/exchange news energies
617 !  Retrieve energy of non-evolving particle(s)
618    if (history_level== 0) then
619      call energies_init(energies)
620      call energies_init(electronpositron%energies_ep)
621      if (dtset%positron>0) energies%e0_electronpositron=etotal_read
622      if (dtset%positron<0) energies%e0_electronpositron=zero
623    else if (history_level== 1) then
624      call energies_init(electronpositron%energies_ep)
625      if (dtset%positron>0) energies%e0_electronpositron=etotal_read
626    else if (history_level== 2) then
627      call energies_copy(energies,electronpositron%energies_ep)
628      call energies_init(energies)
629      energies%e0_electronpositron=electronpositron%e0
630    else if (history_level== 3) then
631      call energies_copy(electronpositron%energies_ep,energies_tmp)
632      call energies_copy(energies,electronpositron%energies_ep)
633      call energies_copy(energies_tmp,energies)
634      energies%e0_electronpositron=electronpositron%e0
635 !    else if (history_level== 4) then
636    end if
637 
638 !  Adjust core psps energy
639    if (electronpositron%calctype/=1) energies%e_corepsp=ecore/ucvol
640 
641 !  -----------------------------------------------------------------------------------------
642 !  Update wavefunctions
643 !  If electronpositron%calctype==1: cg are the positronic WFs, cg_ep are the electronic WFs
644 !  If electronpositron%calctype==2: cg are the electronic WFs, cg_ep are the positronic WFs
645 !  -----------------------------------------------------------------------------------------
646    if (electronpositron%dimcg>0.or.electronpositron%dimcprj>0) then
647 
648      if (history_level==0.or.history_level==1) then
649        electronpositron%cg_ep=zero
650      end if
651 
652      if (history_level==2) then
653        if (electronpositron%dimcg>0) then
654          do icg=1,electronpositron%dimcg
655            electronpositron%cg_ep(1:2,icg)=cg(1:2,icg)
656          end do
657        end if
658        if (dtset%usepaw==1.and.electronpositron%dimcprj>0) then
659          call pawcprj_copy(cprj,electronpositron%cprj_ep)
660        end if
661      end if
662 
663      if (history_level==3) then
664        if (electronpositron%dimcg>0) then
665          do icg=1,electronpositron%dimcg
666            cgtmp(1:2)=electronpositron%cg_ep(1:2,icg)
667            electronpositron%cg_ep(1:2,icg)=cg(1:2,icg)
668            cg(1:2,icg)=cgtmp(1:2)
669          end do
670        end if
671        if (dtset%usepaw==1.and.electronpositron%dimcprj>0) then
672          ABI_ALLOCATE(nlmn,(dtset%natom))
673          ABI_DATATYPE_ALLOCATE(cprj_tmp,(dtset%natom,electronpositron%dimcprj))
674          do iatom=1,dtset%natom
675            nlmn(iatom)=cprj(iatom,1)%nlmn
676          end do
677          call pawcprj_alloc(cprj_tmp,cprj(1,1)%ncpgr,nlmn)
678          ABI_DEALLOCATE(nlmn)
679          call pawcprj_copy(electronpositron%cprj_ep,cprj_tmp)
680          call pawcprj_copy(cprj,electronpositron%cprj_ep)
681          call pawcprj_copy(cprj_tmp,cprj)
682          call pawcprj_free(cprj_tmp)
683          ABI_DATATYPE_DEALLOCATE(cprj_tmp)
684        end if
685      end if
686 
687    end if ! dimcg>0 or dimcprj>0
688 
689 !  -----------------------------------------------------------------------------------------------------------
690 !  Initialize/Update occupations
691 !  If electronpositron%calctype==1: occ are the positronic occupations, occ_ep are the electronic occupations
692 !  If electronpositron%calctype==2: occ are the electronic occupations, occ_ep are the positronic occupations
693 !  -----------------------------------------------------------------------------------------------------------
694 !  When needed, precompute electronic occupations with semiconductor occupancies
695    need_scocc=.false.
696    if (electronpositron%dimocc>0.and.electronpositron%calctype==1.and. &
697 &   (history_level==0.or.history_level==1)) need_scocc=.true.
698    if (electronpositron%calctype==2.and.ireadwf==0.and. &
699 &   (history_level==0.or.history_level==2.or. &
700 &   (history_level==3.and.electronpositron%dimocc==0))) need_scocc=.true.
701    if (need_scocc) then
702      nelect=-dtset%charge
703      do iatom=1,dtset%natom
704        nelect=nelect+dtset%ziontypat(dtset%typat(iatom))
705      end do
706      maxocc=two/real(dtset%nsppol*dtset%nspinor,dp)
707      nocc=(nelect-tol8)/maxocc + 1
708      nocc=min(nocc,dtset%nband(1)*dtset%nsppol)
709      occlast=nelect-maxocc*(nocc-1)
710      ABI_ALLOCATE(scocc,(dtset%nband(1)*dtset%nsppol))
711      scocc=zero
712      if (1<nocc)  scocc(1:nocc-1)=maxocc
713      if (1<=nocc) scocc(nocc)=occlast
714    end if
715 
716 !  ===== PREVIOUS OCCUPATIONS OCC_EP:
717    if (electronpositron%dimocc>0) then
718      if (history_level==0.or.history_level==1) then
719 !      ----- Electronic from scratch
720        if (electronpositron%calctype==1) then
721 !        Initialize electronic occupations with semiconductor occupancies
722          do ikpt=1,dtset%nkpt
723            do iband=1,dtset%nband(1)
724              do isppol=1,dtset%nsppol
725                electronpositron%occ_ep(iband+dtset%nband(1)*(ikpt-1+dtset%nkpt*(isppol-1)))=&
726 &               scocc(isppol+dtset%nsppol*(iband-1))
727              end do
728            end do
729          end do
730        end if
731 !      ----- Positronic from scratch
732        if (electronpositron%calctype==1) then
733 !        Initialize positronic occupations with only one positron (or less)
734          electronpositron%occ_ep(:)=zero
735          isppol=1;iocc=1
736          do ikpt=1,dtset%nkpt
737            electronpositron%occ_ep(iocc)=electronpositron%posocc
738            iocc=iocc+dtset%nband(ikpt+dtset%nkpt*(isppol-1))
739          end do
740        end if
741      end if
742 !    ----- Deduced from occ in memory
743      if (history_level==2) then
744        electronpositron%occ_ep(:)=occ(:)
745      end if
746    end if ! dimocc>0
747 
748 !  ===== CURRENT OCCUPATIONS OCC:
749    if (history_level==0.or.history_level==2.or.(history_level==3.and.electronpositron%dimocc==0)) then
750      if (ireadwf==0) then
751 !      ----- Positronic from scratch
752        if (electronpositron%calctype==1) then
753 !        Initialize positronic occupations with only one positron (or less)
754          occ(:)=zero
755          isppol=1;iocc=1
756          do ikpt=1,dtset%nkpt
757            occ(iocc)=electronpositron%posocc
758            iocc=iocc+dtset%nband(ikpt+dtset%nkpt*(isppol-1))
759          end do
760        end if
761 !      ----- Electronic from scratch
762        if (electronpositron%calctype==2) then
763 !        Initialize electronic occupations with semiconductor occupancies
764          do ikpt=1,dtset%nkpt
765            do iband=1,dtset%nband(1)
766              do isppol=1,dtset%nsppol
767                occ(iband+dtset%nband(1)*(ikpt-1+dtset%nkpt*(isppol-1)))=&
768 &               scocc(isppol+dtset%nsppol*(iband-1))
769              end do
770            end do
771          end do
772        end if
773      end if
774    end if
775 
776 !  ===== EXCHANGE POSITRONIC AND ELECTRONIC OCCUPATIONS (CURRENT AND PREVIOUS)
777    if (history_level==3.and.electronpositron%dimocc>0) then
778      do iocc=1,electronpositron%dimocc
779        occtmp=occ(iocc)
780        occ(iocc)=electronpositron%occ_ep(iocc)
781        electronpositron%occ_ep(iocc)=occtmp
782      end do
783    end if
784 
785    if (need_scocc)  then
786      ABI_DEALLOCATE(scocc)
787    end if
788 
789 !  -----------------------------------------------------------------------------------------------------------
790 !  Initialize/Update eigen energies
791 !  If electronpositron%calctype==1: eigen are the positronic eigen E, eigen_ep are the electronic eigen E
792 !  If electronpositron%calctype==2: eigen are the electronic eigen E, eigen_ep are the positronic eigen E
793 !  -----------------------------------------------------------------------------------------------------------
794 
795 !  ===== PREVIOUS EIGEN ENERGIES EIGEN_EP:
796    if (electronpositron%dimeigen>0) then
797      if (history_level==0.or.history_level==1) then
798 !      ----- Electronic or positronic from scratch
799        electronpositron%eigen_ep(:)=zero
800      end if
801 !    ----- Deduced from eigen in memory
802      if (history_level==2) then
803        electronpositron%eigen_ep(:)=eigen(:)
804      end if
805    end if ! dimeigen>0
806 
807 !  ===== CURRENT EIGEN ENERGIES EIGEN:
808    if (history_level==0.or.history_level==2.or.(history_level==3.and.electronpositron%dimeigen==0)) then
809      if (ireadwf==0) then
810 !      ----- Electronic or positronic from scratch
811        eigen(:)=zero
812      end if
813    end if
814 
815 !  ===== EXCHANGE POSITRONIC AND ELECTRONIC EIGEN ENERGIES (CURRENT AND PREVIOUS)
816    if (history_level==3.and.electronpositron%dimeigen>0) then
817      do iocc=1,electronpositron%dimeigen
818        eigtmp=eigen(iocc)
819        eigen(iocc)=electronpositron%eigen_ep(iocc)
820        electronpositron%eigen_ep(iocc)=eigtmp
821      end do
822    end if
823 
824 !  =============================================
825  end if  ! the type of e-p calculation changes
826 !=============================================
827 
828 !------------------------------------------------------------------
829 !Write messages
830 !------------------------------------------------------------------
831  if (istep_mix==1.and.dtset%positron/=0) then
832 !  Log message
833    if (electronpositron%calctype==0) then
834      message = 'Were are now performing an electronic ground-state calculation...'
835    else if (electronpositron%calctype==1) then
836      message = 'Were are now performing a positronic ground-state calculation...'
837    else if (electronpositron%calctype==2) then
838      message = 'Were are now performing an electronic ground-state calculation in presence of a positron...'
839    end if
840    MSG_COMMENT(message)
841 !  Output message
842    if (dtset%positron<0) then
843      if (electronpositron%calctype==0) then
844        TypeCalcStrg='ELECTRONIC GROUND-STATE CALCULATION'
845      else if (electronpositron%calctype==1) then
846        TypeCalcStrg='POSITRONIC GROUND-STATE CALCULATION IN PRESENCE OF ELECTRONS AND IONS'
847      else if (electronpositron%calctype==2) then
848        TypeCalcStrg='ELECTRONIC GROUND-STATE CALCULATION IN PRESENCE OF A POSITRON'
849      end if
850      if (istep>1) then
851        write(message,'(2a,i3,2a)') ch10,'TC-DFT STEP ',electronpositron%istep,' - ',trim(TypeCalcStrg)
852      else
853        write(message,'(a,i3,2a)') 'TC-DFT STEP ',electronpositron%istep,' - ',trim(TypeCalcStrg)
854      end if
855      call wrtout(ab_out,message,'COLL')
856    end if
857  end if
858 
859 end subroutine setup_positron