./56_xc/rhohxcpositron.F90

ABINIT/rhohxcpositron

ABINIT/rhohxcpositron [ Functions ]

NAME

 rhohxcpositron

FUNCTION

 Calculate the electrons/positron correlation term for the positron

 NOTE

 Copyright (C) 1998-2018 ABINIT group (GJ,MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  gprimd(3,3)=dimensional reciprocal space primitive translations
  mpi_enreg=information about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nhat(nfft,nspden*usepaw)= -PAW only- compensation density
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  nspden=number of spin density components
  n3xccc=dimension of the xccc3d array (0 or nfft).
  paral_kgb=flag for (k,band,FFT) parallelism
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
  ucvol = unit cell volume (Bohr**3)
  usexcnhat= -PAW only- flag controling use of compensation density in Vxc
  usepaw=flag for PAW
  xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
  xc_denpos= lowest allowed density (usually for the computation of the XC functionals)

OUTPUT

  electronpositron%e_xc=electron-positron XC energy
  electronpositron%e_xcdc=Double-counting electron-positron XC energy
  strsxc(6)= contribution of xc to stress tensor (hartree/bohr^3),
  vhartr(nfft)=Hartree potential (returned if option/=0 and option/=10)
  vxcapn=XC electron-positron XC potential for the positron
  vxcavg=unit cell average of Vxc = (1/ucvol) Int [Vxc(r) d^3 r].
  kxcapn(nfft,nkxc)=electron-positron XC kernel (returned only if nkxc/=0)

SIDE EFFECTS

  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation

PARENTS

      energy,rhotov,setvtr

CHILDREN

      mean_fftr,mkdenpos,xcden,xcpositron,xmpi_sum

SOURCE

 55 #if defined HAVE_CONFIG_H
 56 #include "config.h"
 57 #endif
 58 
 59 #include "abi_common.h"
 60 
 61 subroutine rhohxcpositron(electronpositron,gprimd,kxcapn,mpi_enreg,nfft,ngfft,nhat,nkxc,nspden,n3xccc,&
 62 &                         paral_kgb,rhor,strsxc,ucvol,usexcnhat,usepaw,vhartr,vxcapn,vxcavg,xccc3d,xc_denpos)
 63 
 64  use defs_basis
 65  use defs_abitypes
 66  use m_cgtools
 67  use m_xmpi
 68  use m_errors
 69  use m_profiling_abi
 70 
 71  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
 72 
 73 !This section has been created automatically by the script Abilint (TD).
 74 !Do not modify the following lines by hand.
 75 #undef ABI_FUNC
 76 #define ABI_FUNC 'rhohxcpositron'
 77  use interfaces_41_xc_lowlevel
 78  use interfaces_56_xc, except_this_one => rhohxcpositron
 79 !End of the abilint section
 80 
 81  implicit none
 82 
 83 !Arguments ------------------------------------
 84 !scalars
 85  integer,intent(in) :: nfft,nkxc,nspden,n3xccc,paral_kgb,usexcnhat,usepaw
 86  real(dp),intent(in) :: ucvol,xc_denpos
 87  real(dp),intent(out) :: vxcavg
 88  type(electronpositron_type),pointer :: electronpositron
 89 !arrays
 90  integer,intent(in) :: ngfft(18)
 91  real(dp),intent(in) :: gprimd(3,3)
 92  real(dp),intent(in) :: nhat(nfft,nspden*usepaw),rhor(nfft,nspden),xccc3d(n3xccc)
 93  real(dp),intent(out) :: kxcapn(nfft,nkxc),strsxc(6),vhartr(nfft),vxcapn(nfft,nspden)
 94  type(MPI_type),intent(in) :: mpi_enreg
 95 
 96 !Local variables-------------------------------
 97 !scalars
 98  integer :: cplex,ierr,ifft,ishift,iwarn,iwarnp,nfftot,ngr,ngrad,nspden_ep
 99  real(dp) :: exc,excdc,strdiag
100  character(len=500) :: message
101 !arrays
102  real(dp),parameter :: qphon(3)=(/0._dp,0._dp,0._dp/)
103  real(dp) :: vxcavg_tmp(1)
104  real(dp),allocatable :: fxcapn(:),grho2apn(:),rhoe(:,:,:),rhop(:,:),rhotote(:),vxc_ep(:),vxcgr_ep(:)
105 
106 ! *************************************************************************
107 
108  if (electronpositron_calctype(electronpositron)/=1) then
109    message = 'Only electronpositron%calctype=1 allowed !'
110    MSG_BUG(message)
111  end if
112 
113  if (nkxc>3) then
114    message = 'nkxc>3 (Kxc for GGA) not yet implemented !'
115    MSG_ERROR(message)
116  end if
117 
118 !Hartree potential of the positron is zero
119  vhartr=zero
120 
121 !Some allocations/inits
122  ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
123  ngr=0;if (ngrad==2) ngr=nfft
124  ABI_ALLOCATE(fxcapn,(nfft))
125  ABI_ALLOCATE(grho2apn,(ngr))
126  nspden_ep=1;cplex=1;ishift=0
127  iwarn=0;iwarnp=1
128 
129 !Compute total electronic density
130  ABI_ALLOCATE(rhotote,(nfft))
131  rhotote(:)=electronpositron%rhor_ep(:,1)
132  if (n3xccc>0) rhotote(:)=rhotote(:)+xccc3d(:)
133  if (usepaw==1.and.usexcnhat==0) rhotote(:)=rhotote(:)-electronpositron%nhat_ep(:,1)
134 
135 !Extra total electron/positron densities; compute gradients for GGA
136  ABI_ALLOCATE(rhoe,(nfft,nspden_ep,ngrad**2))
137  ABI_ALLOCATE(rhop,(nfft,nspden_ep))
138  call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,paral_kgb,qphon,rhotote,rhoe)
139  if (ngrad==2) grho2apn(:)=rhoe(:,1,2)**2+rhoe(:,1,3)**2+rhoe(:,1,4)**2
140  rhop(:,1)=rhor(:,1);if (usepaw==1.and.usexcnhat==0) rhop(:,1)=rhop(:,1)-nhat(:,1)
141  ABI_DEALLOCATE(rhotote)
142 
143 !Make the densities positive
144  call mkdenpos(iwarn ,nfft,nspden_ep,1,rhoe(:,1,1),xc_denpos)
145  if (.not.electronpositron%posdensity0_limit) then
146    call mkdenpos(iwarnp,nfft,nspden_ep,1,rhop,xc_denpos)
147  end if
148 
149 !Compute electron-positron Vxc_pos, Vxc_el, Fxc, Kxc, ...
150  ABI_ALLOCATE(vxc_ep,(nfft))
151  ABI_ALLOCATE(vxcgr_ep,(ngr))
152  if (nkxc==0) then
153    call xcpositron(fxcapn,grho2apn,electronpositron%ixcpositron,ngr,nfft,electronpositron%posdensity0_limit,&
154 &   rhoe(:,1,1),rhop(:,1),vxc_ep,vxcgr_ep,vxcapn)
155  else
156    call xcpositron(fxcapn,grho2apn,electronpositron%ixcpositron,ngr,nfft,electronpositron%posdensity0_limit,&
157 &   rhoe(:,1,1),rhop(:,1),vxc_ep,vxcgr_ep,vxcapn,dvxce=kxcapn)
158  end if
159  ABI_DEALLOCATE(rhoe)
160  ABI_DEALLOCATE(vxc_ep)
161  ABI_DEALLOCATE(vxcgr_ep)
162  ABI_DEALLOCATE(grho2apn)
163 
164 !Store Vxc and Kxc according to spin components
165  if (nspden>=2) vxcapn(:,2)=vxcapn(:,1)
166  if (nspden==4) vxcapn(:,3:4)=zero
167  if (nkxc==3) then
168    kxcapn(:,1)=two*kxcapn(:,1)
169    kxcapn(:,2)=kxcapn(:,1)
170    kxcapn(:,3)=kxcapn(:,1)
171  end if
172 
173 !Compute XC energies and contribution to stress tensor
174  electronpositron%e_xc  =zero
175  electronpositron%e_xcdc=zero
176  strdiag=zero
177  nfftot=PRODUCT(ngfft(1:3))
178  do ifft=1,nfft
179    electronpositron%e_xc  =electronpositron%e_xc  +fxcapn(ifft)
180    electronpositron%e_xcdc=electronpositron%e_xcdc+vxcapn(ifft,1)*rhor(ifft,1)
181 !  strdiag=strdiag+fxcapn(ifft)   ! Already stored in rhotoxc !
182    strdiag=strdiag-vxcapn(ifft,1)*rhop(ifft,1)
183  end do
184  if (usepaw==1.and.usexcnhat==0) then
185    do ifft=1,nfft
186      electronpositron%e_xcdc=electronpositron%e_xcdc-vxcapn(ifft,1)*nhat(ifft,1)
187    end do
188  end if
189  electronpositron%e_xc  =electronpositron%e_xc  *ucvol/dble(nfftot)
190  electronpositron%e_xcdc=electronpositron%e_xcdc*ucvol/dble(nfftot)
191  strdiag=strdiag/dble(nfftot)
192  ABI_DEALLOCATE(fxcapn)
193  ABI_DEALLOCATE(rhop)
194 
195 !Reduction in case of parallelism
196  if(mpi_enreg%paral_kgb==1)then
197    if(paral_kgb/=0)then
198      exc=electronpositron%e_xc;excdc=electronpositron%e_xcdc
199      call xmpi_sum(exc  ,mpi_enreg%comm_fft,ierr)
200      call xmpi_sum(excdc,mpi_enreg%comm_fft,ierr)
201      electronpositron%e_xc=exc;electronpositron%e_xcdc=excdc
202      call xmpi_sum(strsxc,mpi_enreg%comm_fft,ierr)
203    end if
204  end if
205 
206 !Store stress tensor
207  strsxc(1:3)=strdiag
208  strsxc(4:6)=zero
209 
210 !Compute vxcavg
211  call mean_fftr(vxcapn(:,1),vxcavg_tmp,nfft,nfftot,1)
212  vxcavg=vxcavg_tmp(1)
213 
214 end subroutine rhohxcpositron

TABLE OF CONTENTS

ABINIT/rhohxcpositron [ Functions ]