TABLE OF CONTENTS


ABINIT/xcpositron [ Functions ]

[ Top ] [ Functions ]

NAME

 xcpositron

FUNCTION

 Compute electron-positron correlation potentials and energy density.
 Used electron-positron correlation functional is controlled by ixcpositron argument.
 Returns Fxc, Vxc_pos, Vxc_el from input rhor_pos and rhor_el for positron and electrons.

COPYRIGHT

 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

  grhoe2(ngr)=square of the gradient of electronic density rhoe (needed for GGA)
  ixcpositron=type of electron-positron correlation functional:
     1:  LDA zero positron density limit parametrized by Arponen & Pajanne
         and provided by Boronski & Nieminen [1,2]
     11: LDA zero positron density limit parametrized by Arponen & Pajanne
         and fitted by Sterne & Kaiser [1,3]
     2:  LDA electron-positron correlation
         provided by Puska, Seitsonen, and Nieminen [1,4]
     3:  GGA zero positron density limit parametrized by Arponen & Pajanne
         and provided by Boronski & Nieminen [1,2,5]
     31: GGA zero positron density limit parametrized by Arponen & Pajanne
         and fitted by Sterne & Kaiser [1,3,5]
     See references below
  ngr=size of grho2 array (0 if LDA, npt if GGA)
  npt=number of real space points on which density is provided
  posdensity0_limit=True if we are in the zero positron density limit
  rhoer(npt)=electron density (bohr^-3)
  rhopr(npt)=positron density (bohr^-3)

OUTPUT

  fnxc(npt)=correlation energy per unit volume fxc
  vxce(npt)=correlation potential for electron dfxc/drhoe (hartree)
  vxcp(npt)=correlation potential for positron dfxc/drhop (hartree)
  vxcegr(ngr)= 1/|gradRhoe| dfxc/d|gradRhoe| (empty if LDA, i.e. ngr=0)
  Optional outputs:
    dvxce(npt)=partial second derivatives of the xc energy wr to the electronic density
               dvxce(:)=dVxce/dRhoe
    dvxcp(npt)=partial second derivatives of the xc energy wr to the positronic density
               dvxcp(:)=dVxcp/drhop

NOTES

   References for electron-positron correlation functionals:
         [1] J. Arponen and E. Pajanne, Ann. Phys. (N.Y.) 121, 343 (1979).
         [2] Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986).
         [3] P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991).
         [4] M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994).
         [5] B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1994)

PARENTS

      m_pawxc,rhohxcpositron,rhotoxc

CHILDREN

      invcb

SOURCE

 65 #if defined HAVE_CONFIG_H
 66 #include "config.h"
 67 #endif
 68 
 69 #include "abi_common.h"
 70 
 71 subroutine xcpositron(fnxc,grhoe2,ixcpositron,ngr,npt,posdensity0_limit,rhoer,rhopr,vxce,vxcegr,vxcp,&
 72 &                     dvxce,dvxcp) ! optional arguments
 73 
 74  use defs_basis
 75  use m_errors
 76  use m_profiling_abi
 77 
 78 !This section has been created automatically by the script Abilint (TD).
 79 !Do not modify the following lines by hand.
 80 #undef ABI_FUNC
 81 #define ABI_FUNC 'xcpositron'
 82  use interfaces_41_xc_lowlevel, except_this_one => xcpositron
 83 !End of the abilint section
 84 
 85  implicit none
 86 
 87 !Arguments ------------------------------------
 88 !scalars
 89  integer,intent(in) :: ixcpositron,ngr,npt
 90  logical,intent(in) :: posdensity0_limit
 91 !arrays
 92  real(dp),intent(in) :: grhoe2(ngr),rhoer(npt),rhopr(npt)
 93  real(dp),intent(out) :: fnxc(npt),vxce(npt),vxcegr(ngr),vxcp(npt)
 94  real(dp),intent(out),optional :: dvxce(npt),dvxcp(npt)
 95 
 96 !Local variables-------------------------------
 97 !scalars
 98  integer,parameter :: idebug=0
 99  integer :: ipt
100  logical :: gga,need_dvxce,need_dvxcp
101  real(dp),parameter :: alpha_gga=0.22_dp
102  real(dp),parameter :: ap_a1=-1.56_dp,ap_b1=0.051_dp,ap_c1=-0.081_dp,ap_d1=1.14_dp
103  real(dp),parameter :: ap_a2=-0.92305_dp,ap_b2=-0.05459_dp
104  real(dp),parameter :: ap_a3=-0.6298_dp,ap_b3=-13.15111_dp,ap_c3=2.8655_dp
105  real(dp),parameter :: ap_a4=-179856.2768_dp,ap_b4=186.4207_dp,ap_c4=-0.524_dp
106  real(dp),parameter :: ap_psn_limit=0.7_dp
107  real(dp),parameter :: ap_psn_1=0.9_dp*ap_psn_limit,ap_psn_2=1.1_dp*ap_psn_limit
108  real(dp),parameter :: fpi3=third*four_pi
109  real(dp),parameter :: psn_aa=69.7029_dp,psn_ba=-107.4927_dp,psn_bb=141.8458_dp
110  real(dp),parameter :: psn_ca=23.7182_dp,psn_cb=-33.6472_dp ,psn_cc=5.21152_dp
111  real(dp),parameter :: sk_a=-1.56_dp,sk_b=0.1324_dp,sk_c=-4.092_dp,sk_d=51.96_dp,sk_e=0.7207_dp
112  real(dp),parameter :: rsfac=0.6203504908994000_dp
113  real(dp) :: arse,brse,crse,darse,dbrse,dcrse,d2arse,d2brse,d2crse
114  real(dp) :: d2eps,deps,dexc,dexcdg,dexc_p,d2expgga,dexpgga,d2invrs,dinvrs,d2kf,dkf,d2nqtf2,dnqtf2
115  real(dp) :: drse,drsp,d2exc,d2exc_p,d2rse,d2rsp,d2sqr,dsqr
116  real(dp) :: eexp,eps,exc,exc_p,expgga,invf,dinvf,d2invf,invrhoe,invrhop,invrs,invsqr
117  real(dp) :: kf,logrs,nqtf2,opr2,ratio_ap,ratio_psn,rhoe,rhop,rse,rsp,sqr
118  character(len=500) :: msg
119 !arrays
120  real(dp),allocatable :: rsepts(:),rsppts(:)
121 
122 ! *************************************************************************
123 
124  gga=(ngr==npt)
125  need_dvxce=present(dvxce)
126  need_dvxcp=present(dvxcp)
127 
128  if (gga.and.ixcpositron==2) then 
129    msg = 'xcpositron: GGA not yet implemented for ixcpositron=2 !'
130    MSG_ERROR(msg)
131  end if
132  if (posdensity0_limit.and.ixcpositron==2) then
133    msg = 'xcpositron: ixcpositron=2 cannot be treated in the zero positron density limit !'
134    MSG_ERROR(msg)
135  end if
136  if (abs(ixcpositron)/=1.and.ixcpositron/=11.and.ixcpositron/=2.and.ixcpositron/=3.and.ixcpositron/=31) then
137    msg = 'xcpositron: unknown electron-positron correlation !'
138    MSG_ERROR(msg)
139  end if
140 
141 !Compute density radii for rhor_el, rhor_pos
142  ABI_ALLOCATE(rsepts,(npt))
143  call invcb(rhoer(:),rsepts,npt)
144  rsepts(:)=rsfac*rsepts(:)
145  if (ixcpositron==2) then
146    ABI_ALLOCATE(rsppts,(npt))
147    call invcb(rhopr(:),rsppts,npt)
148    rsppts(:)=rsfac*rsppts(:)
149  end if
150 
151 !Loop over grid points
152 !----------------------------------------------------
153  do ipt=1,npt
154 
155    rhoe=rhoer(ipt)
156    rhop=rhopr(ipt)
157    exc=zero;dexc=zero;d2exc=zero;dexcdg=zero
158 
159    rse=rsepts(ipt)
160    invrhoe=one/rhoe
161    drse=-third*rse*invrhoe
162    if (need_dvxce) d2rse= four/nine*rse*invrhoe**2
163 
164 !  Arponen & Pajane parametrization for electron
165    if (ixcpositron/=11.and.ixcpositron/=31) then
166      if (rse<0.302_dp) then
167        invrs=one/rse;invsqr=sqrt(invrs);logrs=log(rse)
168        exc =ap_a1*invsqr+(ap_b1*logrs+ap_c1)*logrs+ap_d1
169        dexc=drse*invrs*(-half*ap_a1*invsqr+two*ap_b1*logrs+ap_c1)
170        if (need_dvxce) d2exc=(d2rse/drse-drse*invrs)*dexc+drse**2*invrs**2*(quarter*ap_a1*invsqr+two*ap_b1)
171      else if (rse>=0.302_dp.and.rse<=0.56_dp) then
172        invrs=one/rse
173        exc =ap_a2+ap_b2*invrs**2
174        dexc=-drse*ap_b2*two*invrs**3
175        if (need_dvxce) d2exc=d2rse/drse*dexc+six*drse**2*ap_b2*invrs**4
176      else if (rse>0.56_dp.and.rse<=8.0_dp) then
177        invrs=one/(rse+2.5_dp)
178        dinvrs=-drse*invrs**2
179 ! jmb : d2rse initialized only if need_dvxce = .True.
180        if (need_dvxce) d2invrs=-d2rse*invrs**2-two*invrs*drse**2
181        exc =ap_a3+ap_b3*invrs**2+ap_c3*invrs
182        dexc=two*ap_b3*invrs*dinvrs+ap_c3*dinvrs
183        if (need_dvxce) d2exc=two*ap_b3*dinvrs**2+(two*ap_b3*invrs+ap_c3)*d2invrs
184      else
185        exc  =ap_a4*rhoe**2+ap_b4*rhoe+ap_c4
186        dexc =two*ap_a4*rhoe+ap_b4
187        if (need_dvxce) d2exc=two*ap_a4
188      end if
189 
190 !    Sterne & Kaiser parametrization for electron
191    else
192      eexp=exp(-(rse+sk_c)**2/sk_d)
193      opr2=(one+rse**2)
194      arse=atan(rse)
195      exc = sk_a/sqrt(arse)+sk_b*eexp+sk_e
196      dexc= -(two*sk_b*eexp*(sk_c+rse)/sk_d + sk_a/(two*opr2*sqrt(arse)**3))*drse
197      if (need_dvxce) d2exc=-(two*sk_b*eexp*(sk_c+rse)/sk_d + sk_a/(two*opr2*arse**1.5_dp))*d2rse &
198 &     +(two*sk_b*eexp*(two*sk_c**2-sk_d+four*sk_c*rse+two*rse**2)/sk_d**2 &
199 &     +sk_a*(three+four*rse*arse)/(four*opr2**2*sqrt(arse)**5))*drse**2
200    end if
201 
202 !  Puska, Seitsonen and Nieminen parametrization for positron
203    if (ixcpositron==2.and.rse>=ap_psn_1) then
204      rsp=rsppts(ipt)
205      invrhop=one/rhop
206      drsp=-third*rsp*invrhop
207      if (need_dvxcp) d2rsp= four/nine*rsp*invrhop**2
208      exc_p=zero;dexc_p=zero;d2exc_p=zero
209      if (rsp<0.302_dp) then
210        invrs=one/rsp;invsqr=sqrt(invrs);logrs=log(rsp)
211        exc_p =ap_a1*invsqr+(ap_b1*logrs+ap_c1)*logrs+ap_d1
212        dexc_p=drsp*invrs*(-half*ap_a1*invsqr+two*ap_b1*logrs+ap_c1)
213        if (need_dvxcp) d2exc_p=(d2rsp/drsp-drsp*invrs)*dexc_p+drsp**2*invrs**2*(quarter*ap_a1*invsqr+two*ap_b1)
214      else if (rsp>=0.302_dp.and.rsp<=0.56_dp) then
215        invrs=one/rsp
216        exc_p =ap_a2+ap_b2*invrs**2
217        dexc_p=-drsp*ap_b2*two*invrs**3
218        if (need_dvxcp) d2exc_p=d2rsp/drsp*dexc_p+six*drsp**2*ap_b2*invrs**4
219      else if (rsp>0.56_dp.and.rsp<=8.0_dp) then
220        invrs=one/(rsp+2.5_dp)
221        dinvrs=-drsp*invrs**2
222 ! jmb : d2rsp initialized only if need_dvxcp = .True.*
223        if (need_dvxcp) d2invrs=-d2rsp*invrs**2-two*invrs*drsp**2
224        exc_p =ap_a3+ap_b3*invrs**2+ap_c3*invrs
225        dexc_p=two*ap_b3*invrs*dinvrs+ap_c3*dinvrs
226        if (need_dvxcp) d2exc_p=two*ap_b3*dinvrs**2+(two*ap_b3*invrs+ap_c3)*d2invrs
227      else
228        exc_p  =ap_a4*rhop**2+ap_b4*rhop+ap_c4
229        dexc_p =two*ap_a4*rhop+ap_b4
230        if (need_dvxcp) d2exc_p=two*ap_a4
231      end if
232    end if
233 
234 !  GGA correction
235    if (gga) then
236      kf=(three*pi*pi*rhoe)**third
237      nqtf2=(rhoe*sqrt(four*kf/pi))**2
238      eps=grhoe2(ipt)/nqtf2
239      if (eps<zero) then 
240        MSG_ERROR('xcpositron: problem, negative GGA espilon !')
241      end if
242      expgga=exp(-alpha_gga*eps*third)
243 
244      dkf=pi*pi/(sqrt(three*pi*pi*rhoe)**third)
245      d2kf=-two*pi*pi*pi*pi*(three*pi*pi*rhoe)**(-5.0_dp/3.0_dp)
246      sqr=sqrt(four*kf/pi)
247      dsqr=(four*dkf/pi)/(two*sqr)
248      d2sqr=two/(pi*sqr*dkf)*(d2kf*sqr-dsqr*dkf)
249      nqtf2=(rhoe*sqr)**two
250      dnqtf2=two*(sqr+rhoe*dsqr)*rhoe*sqr
251      d2nqtf2=two*(rhoe*sqr*(two*dsqr+rhoe*d2sqr) &
252 &     +sqr*(sqr+rhoe*dsqr) &
253 &     +rhoe*(sqr+rhoe*dsqr) ) 
254      deps=-grhoe2(ipt)*dnqtf2/(nqtf2**two)
255      d2eps=-grhoe2(ipt)/(nqtf2*nqtf2*dnqtf2)*(d2nqtf2*nqtf2*nqtf2-two*nqtf2*dnqtf2*dnqtf2)
256      dexpgga=-alpha_gga*third*deps*expgga
257      d2expgga=-alpha_gga*third*(d2eps*expgga+deps*dexpgga)
258      
259      exc   = exc  *expgga
260      dexc=(dexc*expgga+exc*dexpgga)
261      if (need_dvxce) d2exc=d2exc*expgga+two*dexc*dexpgga+exc*d2expgga
262      if (abs(grhoe2(ipt))<1.e24_dp) dexcdg=-exc*alpha_gga*two_thirds/nqtf2
263    end if
264 
265 !  Computation of XC energy, potentials and kernels
266 !  Zero positron density limit
267    if (ixcpositron/=2.or.rse<ap_psn_1) then
268      fnxc(ipt)=rhop*exc
269      vxce(ipt)=rhop*dexc
270      vxcp(ipt)=exc
271      if (need_dvxce) dvxce(ipt)=rhop*d2exc
272      if (need_dvxcp) dvxcp(ipt)=zero
273      if (gga)       vxcegr(ipt)=rhop*dexcdg
274    else
275 !    Puska, Seitsonen and Nieminen functional
276      arse=psn_aa+psn_ba*rse+psn_ca*rse**2
277      brse=psn_ba+psn_bb*rse+psn_cb*rse**2
278      crse=psn_ca+psn_cb*rse+psn_cc*rse**2
279      darse=(psn_ba+two*psn_ca*rse)*drse
280      dbrse=(psn_bb+two*psn_cb*rse)*drse
281      dcrse=(psn_cb+two*psn_cc*rse)*drse
282      invf=arse+brse*rsp+crse*rsp**2+invrhop/exc+invrhoe/exc_p
283      fnxc(ipt)=one/invf
284      dinvf=darse+dbrse*rsp+dcrse*rsp**2-invrhop*dexc/exc**2-invrhoe**2/exc_p
285      vxce(ipt)=-dinvf/invf**2
286      if (need_dvxce) then
287        d2arse=darse*d2rse/drse+two*psn_ca*drse**2
288        d2brse=dbrse*d2rse/drse+two*psn_cb*drse**2
289        d2crse=dcrse*d2rse/drse+two*psn_cc*drse**2
290        d2invf=d2arse+d2brse*rsp+d2crse*rsp**2 &
291 &       +invrhop*(two*dexc**2/exc-d2exc)/exc**2+two*invrhoe**3/exc_p
292        dvxce(ipt)=(two*dinvf**2/invf-d2invf)/invf**2
293      end if
294      dinvf=(brse+two*crse*rsp)*drsp-invrhop**2/exc-invrhoe*dexc_p/exc_p**2
295      vxcp(ipt)=-dinvf/invf**2
296      if (need_dvxcp) then
297        d2invf=two*crse*drsp+(brse+two*crse*rsp)*d2rsp &
298 &       +two*invrhop**3/exc+invrhoe*(two*dexc_p**2/exc_p-d2exc_p)/exc_p**2
299        dvxcp(ipt)=(two*dinvf**2/invf-d2invf)/invf**2
300      end if
301 !    For small rse, use pure Arponen/Pajanne functional
302 !    Around the limit (rse=0.7, see PSN paper), switch smoothly from PSN to AP
303      if (rse>=ap_psn_1.and.rse<=ap_psn_2) then
304        ratio_psn=(rse-ap_psn_1)/(ap_psn_2-ap_psn_1);ratio_ap=one-ratio_psn
305        fnxc(ipt)=ratio_psn*fnxc(ipt)+ratio_ap*rhop*exc
306        vxce(ipt)=ratio_psn*vxce(ipt)+ratio_ap*rhop*dexc
307        vxcp(ipt)=ratio_psn*vxcp(ipt)+ratio_ap*exc
308        if (need_dvxce) dvxce(ipt)=ratio_psn*dvxce(ipt)+ratio_ap*rhop*d2exc
309        if (need_dvxcp) dvxcp(ipt)=ratio_psn*dvxcp(ipt)
310      end if
311    end if
312 
313 !  Debug statements: use polynomial functionals
314    if (idebug>0) then
315      if (idebug==4) then ! order 4
316        fnxc(ipt)=tol3*((rhop**4+rhoe**4)/12._dp+(rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe**2)
317        vxce(ipt)=tol3*((rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe+rhop*rhoe**2)
318        vxcp(ipt)=tol3*((rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe+rhop*rhoe**2)
319        if (need_dvxce) dvxce(ipt)=tol3*(rhop**3/3._dp+rhop*rhoe**2+rhop**2+two*rhop*rhoe)
320        if (need_dvxcp) dvxcp(ipt)=tol3*(rhoe**3/3._dp+rhoe*rhop**2+rhoe**2+two*rhop*rhoe)
321      end if
322      if (idebug==3) then ! order 3
323        fnxc(ipt)=tol3*((rhop**3+rhoe**3)*third+rhop**2*rhoe+rhop*rhoe**2)
324        vxce(ipt)=tol3*(rhop+rhoe)**2
325        vxcp(ipt)=tol3*(rhop+rhoe)**2
326        if (need_dvxce) dvxce(ipt)=tol3*two*rhoe
327        if (need_dvxcp) dvxcp(ipt)=tol3*two*rhop
328      end if
329      if (idebug==2) then ! order 2
330        fnxc(ipt)=tol3*(rhop+rhoe)**2
331        vxce(ipt)=tol3*two*(rhop+rhoe)
332        vxcp(ipt)=tol3*two*(rhop+rhoe)
333        if (need_dvxce) dvxce(ipt)=tol3*two
334        if (need_dvxcp) dvxcp(ipt)=tol3*two
335      end if
336      if (idebug==1) then ! order 1
337        fnxc(ipt)=tol3*(rhop+rhoe)
338        vxce(ipt)=tol3
339        vxcp(ipt)=tol3
340        if (need_dvxce) dvxce(ipt)=zero
341        if (need_dvxcp) dvxcp(ipt)=zero
342      end if
343    end if
344 
345  end do ! ipt
346 
347  ABI_DEALLOCATE(rsepts)
348  if (ixcpositron==2) then
349    ABI_DEALLOCATE(rsppts)
350  end if
351 
352 !Convert everything in Hartree units
353  fnxc(:)=half*fnxc(:)
354  vxce(:)=half*vxce(:)
355  vxcp(:)=half*vxcp(:)
356  if (need_dvxce) dvxce(:)=half*dvxce(:)
357  if (need_dvxcp) dvxcp(:)=half*dvxcp(:)
358  if (gga)       vxcegr(:)=half*vxcegr(:)
359 
360 end subroutine xcpositron